Drug Discovery in Low Data Regimes: Leveraging a Computational Pipeline for the Discovery of Novel SARS-CoV-2 Nsp14-MTase Inhibitors.

The COVID-19 pandemic, caused by the SARS-CoV-2 virus, has led to significant global morbidity and mortality. A crucial viral protein, the non-structural protein 14 (nsp14), catalyzes the methylation of viral RNA and plays a critical role in viral genome replication and transcription. Due to the low mutation rate in the nsp region among various SARS-CoV-2 variants, nsp14 has emerged as a promising therapeutic target. However, discovering potential inhibitors remains a challenge. In this work, we introduce a computational pipeline for the rapid and efficient identification of potential nsp14 inhibitors by leveraging virtual screening and the NCI open compound collection, which contains 250,000 freely available molecules for researchers worldwide. The introduced pipeline provides a cost-effective and efficient approach for early-stage drug discovery by allowing researchers to evaluate promising molecules without incurring synthesis expenses. Our pipeline successfully identified seven promising candidates after experimentally validating only 40 compounds. Notably, we discovered NSC620333, a compound that exhibits a strong binding affinity to nsp14 with a dissociation constant of 427 ± 84 nM. In addition, we gained new insights into the structure and function of this protein through molecular dynamics simulations. We identified new conformational states of the protein and determined that residues Phe367, Tyr368, and Gln354 within the binding pocket serve as stabilizing residues for novel ligand interactions. We also found that metal coordination complexes are crucial for the overall function of the binding pocket. Lastly, we present the solved crystal structure of the nsp14-MTase complexed with SS148 (PDB:8BWU), a potent inhibitor of methyltransferase activity at the nanomolar level (IC50 value of 70 ± 6 nM). Our computational pipeline accurately predicted the binding pose of SS148, demonstrating its effectiveness and potential in accelerating drug discovery efforts against SARS-CoV-2 and other emerging viruses.

Computational completion of the Aurora interaction region of N-Myc in the Aurora a kinase complex.

Inhibiting protein-protein interactions of the Myc family is a viable pharmacological strategy for modulation of the levels of Myc oncoproteins in cancer. Aurora A kinase (AurA) and N-Myc interaction is one of the most attractive targets of this strategy because formation of this complex blocks proteasomal degradation of N-Myc in neuroblastoma. Two crystallization studies have captured this complex (PDB IDs: 5g1x, 7ztl), partially resolving the AurA interaction region (AIR) of N-Myc. Prompted by the missing N-Myc fragment in these crystal structures, we modeled the complete structure between AurA and N-Myc, and comprehensively analyzed how the incomplete and complete N-Myc behave in complex by molecular dynamics simulations. Molecular dynamics simulations of the incomplete PDB complex (5g1x) repeatedly showed partial dissociation of the short N-Myc fragment (61-89) from the kinase. The missing N-Myc (19-60) fragment was modeled utilizing the N-terminal lobe of AurA as the protein-protein interaction surface, wherein TPX2, a well-known partner of AurA, also binds. Binding free energy calculations along with flexibility analysis confirmed that the complete AIR of N-Myc stabilizes the complex, accentuating the N-terminal lobe of AurA as a binding site for the missing N-Myc fragment (19-60). We further generated additional models consisting of only the missing N-Myc (19-60), and the fused form of TPX2 (7-43) and N-Myc (61-89). These partners also formed more stable interactions with the N-terminal lobe of AurA than did the incomplete N-Myc fragment (61-89) in the 5g1x complex. Altogether, this study provides structural insights into the involvement of the N-terminus of the AIR of N-Myc and the N-terminal lobe of AurA in formation of a stable complex, reflecting its potential for effective targeting of N-Myc.

Computed Protein-Protein Enthalpy Signatures as a Tool for Identifying Conformation Sampling Problems.

Understanding the thermodynamic signature of protein-peptide binding events is a major challenge in computational chemistry. The complexity generated by both components possessing many degrees of freedom poses a significant issue for methods that attempt to directly compute the enthalpic contribution to binding. Indeed, the prevailing assumption has been that the errors associated with such approaches would be too large for them to be meaningful. Nevertheless, we currently have no indication of how well the present methods would perform in terms of predicting the enthalpy of binding for protein-peptide complexes. To that end, we carefully assembled and curated a set of 11 protein-peptide complexes where there is structural and isothermal titration calorimetry data available and then computed the absolute enthalpy of binding. The initial "out of the box" calculations were, as expected, very modest in terms of agreement with the experiment. However, careful inspection of the outliers allows for the identification of key sampling problems such as distinct conformations of peptide termini not being sampled or suboptimal cofactor parameters. Additional simulations guided by these aspects can lead to a respectable correlation with isothermal titration calorimetry (ITC) experiments (R2 of 0.88 and an RMSE of 1.48 kcal/mol overall). Although one cannot know prospectively whether computed ITC values will be correct or not, this work shows that if experimental ITC data are available, then this in conjunction with computed ITC, can be used as a tool to know if the ensemble being simulated is representative of the true ensemble or not. That is important for allowing the correct interpretation of the detailed dynamics of the system with respect to the measured enthalpy. The results also suggest that computational calorimetry is becoming increasingly feasible. We provide the data set as a resource for the community, which could be used as a benchmark to help further progress in this area.

The role of loop dynamics in the prediction of ligand-protein binding enthalpy.

The enthalpic and entropic components of ligand-protein binding free energy reflect the interactions and dynamics between ligand and protein. Despite decades of study, our understanding and hence our ability to predict these individual components remains poor. In recent years, there has been substantial effort and success in the prediction of relative and absolute binding free energies, but the prediction of the enthalpic (and entropic) contributions in biomolecular systems remains challenging. Indeed, it is not even clear what kind of performance in terms of accuracy could currently be obtained for such systems. It is, however, relatively straight-forward to compute the enthalpy of binding. We thus evaluated the performance of absolute enthalpy of binding calculations using molecular dynamics simulation for ten inhibitors against a member of the bromodomain family, BRD4-1, against isothermal titration calorimetry data. Initial calculations, with the AMBER force-field showed good agreement with experiment (R2 = 0.60) and surprisingly good accuracy with an average of root-mean-square error (RMSE) = 2.49 kcal mol-1. Of the ten predictions, three were obvious outliers that were all over-predicted compared to experiment. Analysis of various simulation factors, including parameterization, buffer concentration and conformational dynamics, revealed that the behaviour of a loop (the ZA loop on the periphery of the binding site) strongly dictates the enthalpic prediction. Consistent with previous observations, the loop exists in two distinct conformational states and by considering one or the other or both states, the prediction for the three outliers can be improved dramatically to the point where the R2 = 0.95 and the accuracy in terms of RMSE improves to 0.90 kcal mol-1. However, performance across force-fields is not consistent: if OPLS and CHARMM are used, different outliers are observed and the correlation with the ZA loop behaviour is not recapitulated, likely reflecting parameterization as a confounding problem. The results provide a benchmark standard for future study and comparison.

A recursive molecular docking coupled with energy-based pose-rescoring and MD simulations to identify hsGC ßH-NOX allosteric modulators for cardiovascular dysfunctions.

Modulating the activity of human soluble guanylate cyclase (hsGC) through allosteric regulation of the ßH-NOX domain has been considered as an immediate treatment for cardiovascular disorder (CVDs). Currently available ßH-NOX domain-specific agonists including cinaciguat are unable to deal with the conundrum raised due to oxidative stress in the case of CVDs and their associated comorbidities. Therefore, the idea of investigating novel compounds for allosteric regulation of hsGC activation has been rekindled to circumvent CVDs. Current study aims to identify novel ßH-NOX domain-specific compounds that can selectively turn on sGC functions by modulating the conformational dynamics of the target protein. Through a comprehensive computational drug-discovery approach, we first executed a target-based performance assessment of multiple docking (PLANTS, QVina, LeDock, Vinardo, Smina) scoring functions based on multiple performance metrices. QVina showed the highest capability of selecting true-positive ligands over false positives thus, used to screen 4.8 million ZINC15 compounds against ßH-NOX domain. The docked ligands were further probed in terms of contact footprint and pose reassessment through clustering analysis and PLANTS docking, respectively. Subsequently, energy-based AMBER rescoring of top 100 low-energy complexes, per-residue energy decomposition analysis, and ADME-Tox analysis yielded the top three compounds i.e. ZINC000098973660, ZINC001354120371, and ZINC000096022607. The impact of three selected ligands on the internal structural dynamics of the ßH-NOX domain was also investigated through molecular dynamics simulations. The study revealed potential electrostatic interactions for better conformational dialogue between ßH-NOX domain and allosteric ligands that are critical for the activation of hsGC as compared to the reference compound.

Bridging the Bridging Imidazolate in the Bimetallic Center of the Cu/Zn SOD1 and ALS.

Metallation status of human Cu/Zn superoxide dismutase 1 (SOD1) plays a pivotal role in the pathogenesis of amyotrophic lateral sclerosis (ALS). All of the amino acids found in the bimetallic center have been associated with ALS except for two positions. H63 which forms the bridging imidazolate ion in the bimetallic center and K136 which is not directly involved in coordination but located in the bimetallic center were not reported to be mutated in any of the identified ALS cases. In this study, we investigated the structure and flexibility of five SOD1 variants by using classical molecular dynamics simulations. These variants include three substitutions on the non-ALS-linked positions; H63A, H63R, K136A and ALS-linked positions; G37R, H46R/H48D. We have generated four systems for each variant differing in metallation and presence of the intramolecular disulfide bond. Overall, a total of 24 different dimers including the wild-type were generated and simulated at two temperatures, 298 and 400 K. We have monitored backbone mobility, fluctuations and compactness of the dimer structures to assess whether the hypothetical mutations would behave similar to the ALS-linked variants. Results showed that particularly two mutants, H63R and K136A, drastically affected the dimer dynamics by increasing the fluctuations of the metal binding loops compared with the control mutations. Further, these variants resulted in demetallation of the dimers, highlighting probable ALS toxicity that could be elicited by the SOD1 variants of H63R and K136A. Overall, this study bridges two putative SOD1 positions in the metallic center and ALS, underlining the potential use of atomistic simulations for studying disease variants.

VirtualFlow Ants-Ultra-Large Virtual Screenings with Artificial Intelligence Driven Docking Algorithm Based on Ant Colony Optimization.

The docking program PLANTS, which is based on ant colony optimization (ACO) algorithm, has many advanced features for molecular docking. Among them are multiple scoring functions, the possibility to model explicit displaceable water molecules, and the inclusion of experimental constraints. Here, we add support of PLANTS to VirtualFlow (VirtualFlow Ants), which adds a valuable method for primary virtual screenings and rescoring procedures. Furthermore, we have added support of ligand libraries in the MOL2 format, as well as on the fly conversion of ligand libraries which are in the PDBQT format to the MOL2 format to endow VirtualFlow Ants with an increased flexibility regarding the ligand libraries. The on the fly conversion is carried out with Open Babel and the program SPORES. We applied VirtualFlow Ants to a test system involving KEAP1 on the Google Cloud up to 128,000 CPUs, and the observed scaling behavior is approximately linear. Furthermore, we have adjusted several central docking parameters of PLANTS (such as the speed parameter or the number of ants) and screened 10 million compounds for each of the 10 resulting docking scenarios. We analyzed their docking scores and average docking times, which are key factors in virtual screenings. The possibility of carrying out ultra-large virtual screening with PLANTS via VirtualFlow Ants opens new avenues in computational drug discovery.

Evaluating the Performance of Water Models with Host-Guest Force Fields in Binding Enthalpy Calculations for Cucurbit[7]uril-Guest Systems.

Computational prediction of thermodynamic components with computational methods has become increasingly routine in computer-aided drug design. Although there has been significant recent effort and improvements in the calculation of free energy, the prediction of enthalpy (and entropy) remains underexplored. Furthermore, there has been relatively little work reported so far that attempts to comparatively assess how well different force fields and water models perform in conjunction with each other. Here, we report a comprehensive assessment of force fields and water models using host-guest systems that mimic many features of protein-ligand systems. These systems are computationally inexpensive, possibly because of their small size compared to protein-ligand systems. We present absolute enthalpy calculations using the multibox approach on a set of 25 cucurbit[7]uril-guest pairs. Eight water models were considered (TIP3P, TIP4P, TIP4P-Ew, SPC, SPC/E, OPC, TIP5P, Bind3P), along with five force fields commonly used in the literature (GAFFv1, GAFFv2, CGenFF, Parsley, and SwissParam). We observe that host-guest binding enthalpies are strongly sensitive to the selection of force field and water model. In terms of water models, we find that TIP3P and its derivative Bind3P are the best performing models for this particular host-guest system. The performance is generally better for aliphatic compounds than for aromatic ones, suggesting that aromaticity remains a difficult property to include accurately in these simple force fields.

Exploring the structural basis of conformational heterogeneity and autoinhibition of human cGMP-specific protein kinase Iα through computational modelling and molecular dynamics simulations.

Protein kinase Iα (PKGIα) is a pivotal cyclic guanosine monophosphate (cGMP) signalling protein. Major steps related to the structural plasticity of PKGIα have been inferred but the structural aspects of the auto-inhibition and multidomain tertiary organization of human PKGIα in active and inactive form are not clear. Here we combine computational comparative modelling, protein-protein docking and molecular dynamics (MD) simulations to investigate structural details of the repressed state of the catalytic domain of PKGIα. Exploration of the potential inhibitory conformation of the auto-inhibitory domain (AI) within the catalytic cleft reveals that the pseudo-substrate motif binds with residues of the glycine rich loop and substrate-binding lobe. Dynamic changes as a result of coupling of the catalytic and AI domains are also investigated. The three-dimensional homodimeric models of PKGIα in the active and inactive state indicate that PKGIα in its inactive-state attains a compact globular structure where cyclic nucleotide binding (CNB-A/B) domains are buried, whereas the catalytic domains are inaccessible with their substrate-binding pockets facing the N-terminal of CNB-A. Contrary to this, the active-state model of PKGIα shows an extended conformation where CNB-A/B domains are slightly rearranged and the catalytic domains of homodimer flanking the C-terminal with their substrate binding lobes free to entrap downstream proteins. These findings are consistent with previously reported static images of the multidomain organization of PKGIα. Structural insights pertaining to the conformational heterogeneity and auto-inhibition of PKGIα provided in this study may help to understand the dynamics-driven effective regulation of PKGIα.

Comprehensive evaluation of the MM-GBSA method on bromodomain-inhibitor sets.

MM-PB/GBSA methods represent a higher-level scoring theory than docking. This study reports an extensive testing of different MM-GBSA scoring schemes on two bromodomain (BRD) datasets. The first set is composed of 24 BRPF1 complexes, and the second one is a nonredundant set constructed from the PDBbind and composed of 28 diverse BRD complexes. A variety of MM-GBSA schemes were analyzed to evaluate the performance of four protocols with different numbers of minimization and MD steps, 10 different force fields and three different water models. Results showed that neither additional MD steps nor unfixing the receptor atoms improved scoring or ranking power. On the contrary, our results underscore the advantage of fixing receptor atoms or limiting the number of MD steps not only for a reduction in the computational costs but also for boosting the prediction accuracy. Among Amber force fields tested, ff14SB and its derivatives rather than ff94 or polarized force fields provided the most accurate scoring and ranking results. The TIP3P water model yielded the highest scoring and ranking power compared to the others. Posing power was further evaluated for the BRPF1 set. A slightly better posing power for the protocol which uses both minimization and MD steps with a fixed receptor than the one which uses only minimization with a fully flexible receptor-ligand system was observed. Overall, this study provides insights into the usage of the MM-GBSA methods for screening of BRD inhibitors, substantiating the benefits of shorter protocols and latest force fields and maintaining the crystal waters for accuracy.

Comparative Assessment of Seven Docking Programs on a Nonredundant Metalloprotein Subset of the PDBbind Refined.

Extensive usage of molecular docking for computer-aided drug discovery resulted in development of numerous programs with versatile scoring and posing algorithms. Selection of the docking program among these vast number of options is central to the outcome of drug discovery. To this end, comparative assessment studies of docking offer valuable insights into the selection of the optimal tool. Despite the availability of various docking assessment studies, the performance difference of docking programs has not been well addressed on metalloproteins which comprise a substantial portion of the human proteome and have been increasingly targeted for treatment of a wide variety of diseases. This study reports comparative assessment of seven docking programs on a diverse metalloprotein set which was compiled for this study. The refined set of the PDBbind (2017) was screened to gather 710 complexes with metal ion(s) closely located to the ligands (<4 Å). The redundancy was eliminated by clustering and overall 213 complexes were compiled as the nonredundant metalloprotein subset of the PDBbind refined. The scoring, ranking, and posing powers of seven noncommercial docking programs, namely, AutoDock4, AutoDock4Zn, AutoDock Vina, Quick Vina 2, LeDock, PLANTS, and UCSF DOCK6, were comprehensively evaluated on this nonredundant set. Results indicated that PLANTS (80%) followed by LeDock (77%), QVina (76%), and Vina (73%) had the most accurate posing algorithms while AutoDock4 (48%) and DOCK6 (56%) were the least successful in posing. Contrary to their moderate-to-high level of posing success, none of the programs was successful in scoring or ranking of the binding affinities (r2 ≈ 0). Screening power was further evaluated by using active-decoy ligand sets for a large compilation of metalloprotein targets. PLANTS stood out among other programs to be able to enrich the active ligand for every target, underscoring its robustness for screening of metalloprotein inhibitors. This study provides useful information for drug discovery studies targeting metalloproteins.

Insights into an alternative benzofuran binding mode and novel scaffolds of polyketide synthase 13 inhibitors.

Small molecules targeting biosynthesis of mycolic acids in the tuberculosis causing bacterium carry high potential for anti-tuberculosis drug discovery. Hitherto, benzofuran containing molecules were identified to target the thioesterase domain of polyketide synthase 13 (Pks13), one of the crucial constituents of this pathway. Among these, TAM16 was also reported to be highly potent in vivo. Here we performed a multi-stage virtual screening methodology recruiting both ligand- and structure-based screening tools for identification of novel Pks13 inhibitors. The large ZINC15 database comprising more than 1 billion ligands was reduced to 21,277 ligands with benzofuran rings and similarity to TAM16. This collection was screened by docking and the 20 top scoring ligands were further analyzed by molecular dynamics simulations. Molecular mechanics/generalized Born surface area (MM-GBSA) based binding free energy predictions confirmed five molecules in the ZINC15 database lead to remarkable increases in the binding free energy of TAM16. Essentially, the most potent ligand ZINC840169872, which carries a distinct scaffold with a cyclohexane group fused to the furan ring, produced a twofold change in the ligand efficiency of TAM16. Further, assessments of the structure-based tools on five different Pks13-TAM complex structures suggested a high-level agreement with the experiments, substantiating the validity of our methodology for screening Pks13 inhibitors. Overall, these in silico insights into a low energy benzofuran-based scaffold and an alternative binding mode for the benzofuran ring unravel viable strategies to generate potent anti-tuberculosis drugs, accentuating the applications of virtual screening approaches for exploring large compound libraries that cannot be easily accessed by experimentation. Graphical abstract In silico virtual screening of ZINC15 database for the identification of potential inhibitors targeting PKS13.

In silico identification of inhibitors targeting N-Terminal domain of human Replication Protein A.

Replication Protein A (RPA) mediates DNA Damage Response (DDR) pathways through protein-protein interactions (PPIs). Targeting the PPIs formed between RPA and other DNA Damage Response (DDR) mediators has become an intriguing area of research for cancer drug discovery. A number of studies applied different methods ranging from high throughput screening approaches to fragment-based drug design tools to discover RPA inhibitors. Although these methods are robust, virtual screening approaches may be allocated as an alternative to such experimental methods, especially for screening of large libraries. Here we report the comprehensive screening of the large database, ZINC15 composed of ∼750 M compounds and the comparison of the identified ligands with the previously known inhibitors by means of binding affinity and drug-likeness. Initially, a ligand library sharing similarity with a promising inhibitor of the N-terminal domain of the RPA70 subunit (RPA70N) was generated by screening of the ZINC15 library. 46,999 ligands were collected and screened by LeDock which produced a satisfactory correlation with the experimental values (R2 = 0.77). 10 of the top-scoring ligands in LeDock were directly progressed to molecular dynamics (MD) simulations, while 10 additional ligands were also selected based on their LeDock scores and the presence of a functional group that could interact with the key amino acids in the RPA70N cleft. MD simulations were used to predict the binding free energy of the ligands by the MM-PBSA method which produced a high level of agreement with the experiments (R2 = 0.85). Binding free energy predictions pointed out 2 ligands with higher binding affinity than any of the reference inhibitors. Particularly the ligand ZINC000753854163 exhibited superior drug-likeness features than any of the known inhibitors. Overall, this study reports ZINC000753854163 as a possible inhibitor of RPA70N, reflecting its possible use in RPA70N targeted cancer therapy.

An integrated approach towards the development of novel antifungal agents containing thiadiazole: synthesis and a combined similarity search, homology modelling, molecular dynamics and molecular docking study.

BACKGROUND: This study aims to synthesise and characterise novel compounds containing 2-amino-1,3,4-thiadiazole and their acyl derivatives and to investigate antifungal activities. Similarity search, molecular dynamics and molecular docking were also studied to find out a potential target and enlighten the inhibition mechanism. RESULTS: As a first step, 2-amino-1,3,4-thiadiazole derivatives (compounds 3 and 4) were synthesised with high yields (81 and 84%). The target compounds (6a-n and 7a-n) were then synthesised with moderate to high yields (56-87%) by reacting 3 and 4 with various acyl chloride derivatives (5a-n). The synthesized compounds were characterized using the IR, 1H-NMR, 13C-NMR, Mass, X-ray (compound 7n) and elemental analysis techniques. Later, the in vitro antifungal activities of the synthesised compounds were determined. The inhibition zones exhibited by the compounds against the tested fungi, their minimum fungicidal activities, minimum inhibitory concentration and the lethal dose values (LD50) were determined. The compounds exhibited moderate to high levels of activity against all tested pathogens. Finally, in silico modelling was used to enlighten inhibition mechanism using ligand and structure-based methods. As an initial step, similarity search was carried out and the resulting proteins that belong to Homo sapiens were used as reference in sequence similarity search to find the corresponding amino acid sequences in target organisms. Homology modelling was used to construct the protein structure. The stabilised protein structure obtained from molecular dynamics simulation was used in molecular docking. CONCLUSION: The overall results presented here might be a good starting point for the identification of novel and more active compounds as antifungal agents.