Quantum simulations of SARS-CoV-2 main protease Mpro enable high-quality scoring of diverse ligands.

The COVID-19 pandemic has led to unprecedented efforts to identify drugs that can reduce its associated morbidity/mortality rate. Computational chemistry approaches hold the potential for triaging potential candidates far more quickly than their experimental counterparts. These methods have been widely used to search for small molecules that can inhibit critical proteins involved in the SARS-CoV-2 replication cycle. An important target is the SARS-CoV-2 main protease Mpro, an enzyme that cleaves the viral polyproteins into individual proteins required for viral replication and transcription. Unfortunately, standard computational screening methods face difficulties in ranking diverse ligands to a receptor due to disparate ligand scaffolds and varying charge states. Here, we describe full density functional quantum mechanical (DFT) simulations of Mpro in complex with various ligands to obtain absolute ligand binding energies. Our calculations are enabled by a new cloud-native parallel DFT implementation running on computational resources from Amazon Web Services (AWS). The results we obtain are promising: the approach is quite capable of scoring a very diverse set of existing drug compounds for their affinities to M pro and suggest the DFT approach is potentially more broadly applicable to repurpose screening against this target. In addition, each DFT simulation required only ~ 1 h (wall clock time) per ligand. The fast turnaround time raises the practical possibility of a broad application of large-scale quantum mechanics in the drug discovery pipeline at stages where ligand diversity is essential.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery.

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Rigorous Free Energy Simulations in Virtual Screening.

Virtual high throughput screening (vHTS) in drug discovery is a powerful approach to identify hits: when applied successfully, it can be much faster and cheaper than experimental high-throughput screening approaches. However, mainstream vHTS tools have significant limitations: ligand-based methods depend on knowledge of existing chemical matter, while structure-based tools such as docking involve significant approximations that limit their accuracy. Recent advances in scientific methods coupled with dramatic speedups in computational processing with GPUs make this an opportune time to consider the role of more rigorous methods that could improve the predictive power of vHTS workflows. In this Perspective, we assert that alchemical binding free energy methods using all-atom molecular dynamics simulations have matured to the point where they can be applied in virtual screening campaigns as a final scoring stage to prioritize the top molecules for experimental testing. Specifically, we propose that alchemical absolute binding free energy (ABFE) calculations offer the most direct and computationally efficient approach within a rigorous statistical thermodynamic framework for computing binding energies of diverse molecules, as is required for virtual screening. ABFE calculations are particularly attractive for drug discovery at this point in time, where the confluence of large-scale genomics data and insights from chemical biology have unveiled a large number of promising disease targets for which no small molecule binders are known, precluding ligand-based approaches, and where traditional docking approaches have foundered to find progressible chemical matter.

Predicting mutations deleterious to function in beta-lactamase TEM1 using MM-GBSA.

Missense mutations can have disastrous effects on the function of a protein. And as a result, they have been implicated in numerous diseases. However, the majority of missense variants only have a nominal impact on protein function. Thus, the ability to distinguish these two classes of missense mutations would greatly aid drug discovery efforts in target identification and validation as well as medical diagnosis. Monitoring the co-occurrence of a given missense mutation and a disease phenotype provides a pathway for classifying functionally disrupting missense mutations. But, the occurrence of a specific missense variant is often extremely rare making statistical links challenging to infer. In this study, we benchmark a physics-based approach for predicting changes in stability, MM-GBSA, and apply it to classifying mutations as functionally disrupting. A large and diverse dataset of 990 residue mutations in beta-lactamase TEM1 is used to assess performance as it is rich in both functionally disrupting mutations and functionally neutral/beneficial mutations. On this dataset, we compare the performance of MM-GBSA to alternative strategies for predicting functionally disrupting mutations. We observe that the MM-GBSA method obtains an area under the curve (AUC) of 0.75 on the entire dataset, outperforming all other predictors tested. More importantly, MM-GBSA's performance is robust to various divisions of the dataset, speaking to the generality of the approach. Though there is one notable exception: Mutations on the surface of the protein are the mutations that are the most difficult to classify as functionally disrupting for all methods tested. This is likely due to the many mechanisms available to surface mutations to disrupt function, and thus provides a direction of focus for future studies.

Physics-based enzyme design: predicting binding affinity and catalytic activity.

Computational enzyme design is an emerging field that has yielded promising success stories, but where numerous challenges remain. Accurate methods to rapidly evaluate possible enzyme design variants could provide significant value when combined with experimental efforts by reducing the number of variants needed to be synthesized and speeding the time to reach the desired endpoint of the design. To that end, extending our computational methods to model the fundamental physical-chemical principles that regulate activity in a protocol that is automated and accessible to a broad population of enzyme design researchers is essential. Here, we apply a physics-based implicit solvent MM-GBSA scoring approach to enzyme design and benchmark the computational predictions against experimentally determined activities. Specifically, we evaluate the ability of MM-GBSA to predict changes in affinity for a steroid binder protein, catalytic turnover for a Kemp eliminase, and catalytic activity for α-Gliadin peptidase variants. Using the enzyme design framework developed here, we accurately rank the most experimentally active enzyme variants, suggesting that this approach could provide enrichment of active variants in real-world enzyme design applications.

Structure-based approach to the prediction of disulfide bonds in proteins.

Protein engineering remains an area of growing importance in pharmaceutical and biotechnology research. Stabilization of a folded protein conformation is a frequent goal in projects that deal with affinity optimization, enzyme design, protein construct design, and reducing the size of functional proteins. Indeed, it can be desirable to assess and improve protein stability in order to avoid liabilities such as aggregation, degradation, and immunogenic response that may arise during development. One way to stabilize a protein is through the introduction of disulfide bonds. Here, we describe a method to predict pairs of protein residues that can be mutated to form a disulfide bond. We combine a physics-based approach that incorporates implicit solvent molecular mechanics with a knowledge-based approach. We first assign relative weights to the terms that comprise our scoring function using a genetic algorithm applied to a set of 75 wild-type structures that each contains a disulfide bond. The method is then tested on a separate set of 13 engineered proteins comprising 15 artificial stabilizing disulfides introduced via site-directed mutagenesis. We find that the native disulfide in the wild-type proteins is scored well, on average (within the top 6% of the reasonable pairs of residues that could form a disulfide bond) while 6 out of the 15 artificial stabilizing disulfides scored within the top 13% of ranked predictions. Overall, this suggests that the physics-based approach presented here can be useful for triaging possible pairs of mutations for disulfide bond formation to improve protein stability.

AAK1 identified as an inhibitor of neuregulin-1/ErbB4-dependent neurotrophic factor signaling using integrative chemical genomics and proteomics.

Target identification remains challenging for the field of chemical biology. We describe an integrative chemical genomic and proteomic approach combining the use of differentially active analogs of small molecule probes with stable isotope labeling by amino acids in cell culture-mediated affinity enrichment, followed by subsequent testing of candidate targets using RNA interference-mediated gene silencing. We applied this approach to characterizing the natural product K252a and its ability to potentiate neuregulin-1 (Nrg1)/ErbB4 (v-erb-a erythroblastic leukemia viral oncogene homolog 4)-dependent neurotrophic factor signaling and neuritogenesis. We show that AAK1 (adaptor-associated kinase 1) is a relevant target of K252a, and that the loss of AAK1 alters ErbB4 trafficking and expression levels, providing evidence for a previously unrecognized role for AAK1 in Nrg1-mediated neurotrophic factor signaling. Similar strategies should lead to the discovery of novel targets for therapeutic development.

FURSMASA: a new approach to rapid scoring functions that uses a MD-averaged potential energy grid and a solvent-accessible surface area term with parameters GA fit to experimental data.

We demonstrate a new approach to the development of scoring functions through the formulation and parameterization of a new function, which can be used both for rapidly ranking the binding of ligands to proteins and for estimating relative aqueous molecular solubilities. The intent of this work is to introduce a new paradigm for creation of scoring functions, wherein we impose the following criteria upon the function: (1) simple; (2) intuitive; (3) requires no postparameterization tweaking; (4) can be applied (without reparameterization) to multiple target systems; and (5) can be rapidly evaluated for any potential ligand. Following these criteria, a new function, FURSMASA (function for rapid scoring using an MD-averaged grid and the accessible surface area) has been developed. Three novel features of the function include: (1) use of an MD-averaged potential energy grid for ligand-protein interactions, rather than a simple static grid; (2) inclusion of a term that depends on the change in the solvent-accessible surface area changes on an atomic (not molecular) basis; and (3) use of the recently derived predictive index (PI) target when optimizing the function, which focuses the function on its intended purpose of relative ranking. A genetic algorithm is used to optimize the function against test data sets that include ligands for the following proteins: IMPDH, p38, gyrase B, HIV-1, and TACE, as well as the Syracuse Research solubility database. We find that the function is predictive, and can simultaneously fit all the test data sets with cross-validated predictive indices ranging from 0.68 to 0.82. As a test of the ability of this function to predict binding for systems not in the training set, the resulting fitted FURSAMA function is then applied to 23 ligands of the COX-2 enzyme. Comparing the results for COX-2 against those obtained using a variety of well-known rapid scoring functions demonstrates that FURSMASA outperforms all of them in terms of the PI and correlation coefficient. We also find that the FURSAMA function is able to reliably predict the water solubility for 1032 compounds from the Syracuse Research solubility database with a cross-correlated PI of 0.84 and a correlation coefficient R(2) of 0.69. This prediction, which is based solely on a term derived from the atom-based solvent-accessible surface areas, compares favorably with the best prediction methods in the literature, most of which are more complex and/or require experimental data. Finally, as a rigorous test of the applicability to database screening, we apply FURSMASA to large active/decoy ligand databases for IMPDH (400 actives vs. 10,000 decoys), p38 (502 actives vs. 10,000 decoys), and HIV (787 actives vs. 10,000 decoys) used in earlier work to critically evaluate many popular scoring functions, and find that FURSMASA performs surprisingly well for IMPDH and HIV.

Evaluating the molecular mechanics poisson-boltzmann surface area free energy method using a congeneric series of ligands to p38 MAP kinase.

The recently described molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method for calculating free energies is applied to a congeneric series of 16 ligands to p38 MAP kinase whose binding constants span approximately 2 orders of magnitude. These compounds have previously been used to test and compare other free energy calculation methods, including thermodynamic integration (TI), OWFEG, ChemScore, PLPScore, and Dock Energy Score. We find that the MM-PBSA performs relatively poorly for this set of ligands, yielding results much inferior to those from TI or OWFEG, inferior to Dock Energy Score, and not appreciably better than ChemScore or PLPScore but at an appreciably larger computational cost than any of these other methods. This suggests that one should be selective in applying the MM-PBSA method and that for systems that are amenable to other free energy approaches, these other approaches may be preferred. We also examine the single simulation approximation for MM-PBSA, whereby the required ligand and protein trajectories are extracted from a single MD simulation rather than two separate MD runs. This assumption, sometimes used to speed the MM-PBSA calculation, is found to yield significantly inferior results with only a moderate net percentage reduction in total simulation time.