Accurately measuring compound binding affinities is key to driving the pharmaceutical development process. Rigorous physics-based in silico approaches, particularly alchemical free energy methods, have become a gold-standard tool for estimating compound affinity changes. Here we present the results of a large-scale precompetitive collaborative assessment of relative binding free energy (RBFE) calculations generated by 15 pharmaceutical companies. We evaluate an open-source and MIT-licensed RBFE protocol from the Open Free Energy (OpenFE) ecosystem across both public and blinded private datasets, encompassing over 1,700 ligands in total. For the public dataset, the weighted RMSE across the 58 systems was 1.73(1.53)(1.96) kcal/mol, with 10 of the systems reaching sub-kcal/mol accuracy. For the private dataset, the weighted RMSE across the 37 systems was 2.44(1.94)(3.06) kcal/mol, with only 2 of the systems reaching sub-kcal/mol accuracy, reflecting the increased complexity of real-world drug discovery. The protocol's performance was system-dependent, with no single dominant error source, indicating that accuracy is primarily influenced by input quality and transformation type. Overall, these benchmark results are encouraging and indicate that OpenFE is ready for large-scale industrial applications, with an "out-of-the-box" accuracy that approaches that of commercial solutions. While comparison against published FEP+ results, which were obtained after manual parameter optimization, shows comparable ranking statistics, a gap remains in error statistics, for which we outline possible paths toward improvement. The protocol meets key criteria required for production use in an industrial setting: it shows robust performance, generates reproducible results, and achieves both sufficient throughput and rapid convergence.

Recent advances in technologies like cryoEM structure resolution and protein de novo folding prediction have resulted in a wealth of macromolecular structures that have not been resolved to the level of detail a high-resolution X-ray crystal structure could provide. Taking full advantage of these structures for rational drug design would benefit from additional validation and refinement. In this presentation, we investigate if computational refinement and structure-based modeling methods can be utilized to generate reliable complex poses. We present a solution to the induced fit docking problem for protein−ligand binding by combining ligand-based pharmacophore docking, rigid receptor docking, and protein structure prediction with explicit solvent molecular dynamics simulations. This methodology succeeded in determining protein−ligand binding modes with a root-mean-square deviation within 2.5 A compared to experiment in over 90% of cross-docking cases in our testing. Applications of the predicted ligand-receptor structure in free energy perturbation calculations for additional validation is demonstrated.

In this study, we generated a matched molecular pair dataset of halogen/deshalogen compounds with reliable binding affinity data and structural binding mode information from public databases. The workflow includes automated system preparation and setup of free energy perturbation relative binding free energy calculations. We demonstrate the suitability of these datasets to investigate the performance of molecular mechanics force fields and molecular simulation algorithms for the purpose of in silico affinity predictions in lead optimization. Our datasets of a total of 115 matched molecular pairs show highly accurate binding free energy predictions with an average error of <1 kcal/mol despite the semi-automated calculation scheme. We quantify the accuracy of the optimized potential for liquid simulations (OPLS) force field to predict the effect of halogen addition to compounds, a commonly employed chemical modification in the design of drug-like molecules.

Here we present an evaluation of the binding affinity prediction accuracy of the free energy calculation method FEP+ on internal active drug discovery projects and on a large new public benchmark set.

Hypophosphatasia (HPP) is a rare metabolic disorder characterized by low tissue-nonspecific alkaline phosphatase (TNSALP) typically caused by ALPL gene mutations. HPP is heterogeneous, with clinical presentation correlating with residual TNSALP activity and/or dominant-negative effects (DNE). We measured residual activity and DNE for 155 ALPL variants by transient transfection and TNSALP enzymatic activity measurement. Ninety variants showed low residual activity and 24 showed DNE. These results encompass all missense variants with carrier frequencies above 1/25,000 from the Genome Aggregation Database. We used resulting data as a reference to develop a new computational algorithm that scores ALPL missense variants and predicts high/low TNSALP enzymatic activity. Our approach measures the effects of amino acid changes on TNSALP dimer stability with a physics-based implicit solvent energy model. We predict mutation deleteriousness with high specificity, achieving a true-positive rate of 0.63 with false-positive rate of 0, with an area under receiver operating curve (AUC) of 0.9, better than all in silico predictors tested. Combining this algorithm with other in silico approaches can further increase performance, reaching an AUC of 0.94. This study expands our understanding of HPP heterogeneity and genotype/phenotype relationships with the aim of improving clinical ALPL variant interpretation.

Optimizing the solubility of small molecules is important in a wide variety of contexts, including in drug discovery where the optimization of aqueous solubility is often crucial to achieve oral bioavailability. In such a context, solubility optimization cannot be successfully pursued by indiscriminate increases in polarity, which would likely reduce permeability and potency. Moreover, increasing polarity may not even improve solubility itself in many cases, if it stabilizes the solid-state form. Here we present a novel physics-based approach to predict the solubility of small molecules, that takes into account three-dimensional solid-state characteristics in addition to polarity. The calculated solubilities are in good agreement with experimental solubilities taken both from the literature as well as from several active pharmaceutical discovery projects. This computational approach enables strategies to optimize solubility by disrupting the three-dimensional solid-state packing of novel chemical matter, illustrated here for an active medicinal chemistry campaign.

Numerous types of quantum chemical calculations and protocols have been successfully applied to computing of small, uncomplicated organic molecules. Here, we argue for the need to shift attention to more challenging molecules that are marked by an interplay of complicating factors such as conformational, tautomeric, steric, and other effects. The challenge is not in choosing the right quantum chemical method and solvation model but in combining the existing methods to simultaneously and accurately describe the breadth of chemical and physical phenomena that give rise to the experimentally observed . The complexity of the phenomena that must be considered begs for the need for a greater automation of prediction workflows. We review our experience with these challenges and outline paths for future progress in the direction of tackling prediction of complex organic molecules.

The membrane alignment of helical amphiphilic peptides in oriented phospholipid bilayers can be obtained as ensemble and time averages from solid state 2H NMR by fitting the quadrupolar splittings to ideal α-helices. At the same time, molecular dynamics (MD) simulations can provide atomistic insight into peptide-membrane systems. Here, we evaluate the potential of MD simulations to complement the experimental NMR data that is available on three exemplary systems: the natural antimicrobial peptide PGLa and the two designer-made peptides MSI-103 and KIA14, whose sequences were derived from PGLa. Each peptide was simulated for 1 μs in a DMPC lipid bilayer. We calculated from the MD simulations the local angles which define the side chain geometry with respect to the peptide helix. The peptide orientation was then calculated (i) directly from the simulation, (ii) from back-calculated MD-derived NMR splittings, and (iii) from experimental 2H NMR splittings. Our findings are that (1) the membrane orientation and secondary structure of the peptides found in the NMR analysis are generally well reproduced by the simulations; (2) the geometry of the side chains with respect to the helix backbone can deviate significantly from the ideal structure depending on the specific residue, but on average all side chains have the same orientation; and (3) for all of our peptides, the azimuthal rotation angle found from the MD-derived splittings is about 15° smaller than the experimental value.

The therapeutic effect of targeted kinase inhibitors can be significantly reduced by intrinsic or acquired resistance mutations that modulate the affinity of the drug for the kinase. In cancer, the majority of missense mutations are rare, making it difficult to predict their impact on inhibitor affinity. This complicates the practice of precision medicine, pairing of patients with clinical trials, and development of next-generation inhibitors. Here, we examine the potential for alchemical free-energy calculations to predict how kinase mutations modulate inhibitor affinities to Abl, a major target in chronic myelogenous leukemia (CML). We find these calculations can achieve useful accuracy in predicting resistance for a set of eight FDA-approved kinase inhibitors across 144 clinically-identified point mutations, achieving a root mean square error in binding free energy changes of 1.10.91.3 kcal/mol (95% confidence interval) and correctly classifying mutations as resistant or susceptible with 888293% accuracy. Since these calculations are fast on modern GPUs, this benchmark establishes the potential for physical modeling to collaboratively support the rapid assessment and anticipation of the potential for patient mutations to affect drug potency in clinical applications.

We present a molecular dynamics simulation study of alkali metal cation transport through the double‐helical and the head‐to‐head conformers of the gramicidin ion channel. Our approach is based on a thermodynamic integration network, which consists of a sequence of transport reactions, absolute free energies of solvation and cycles of alchemical transmutations of the ions. In this manner, we can reliably estimate free energies and their statistical errors via a least‐squares method without imposing external forces on the system. Within the double helical channel, we find a free energy surface typical for hopping transport between isoenergetic sites of ion localization, separated by comparatively large activation barriers. For fast transport through the head‐to‐head conformation, the thermodynamic network scheme starts to break down. © 2018 Wiley Periodicals, Inc.

Optimization of fragment size d-amino acid oxidase (DAAO) inhibitors was investigated using a combination of computational and experimental methods. Retrospective free energy perturbation (FEP) calculations were performed for benzo[d]isoxazole derivatives, a series of known inhibitors with two potential binding modes derived from X-ray structures of other DAAO inhibitors. The good agreement between experimental and computed binding free energies in only one of the hypothesized binding modes strongly support this bioactive conformation. Then, a series of 1-H-indazol-3-ol derivatives formerly not described as DAAO inhibitors was investigated. Binding geometries could be reliably identified by structural similarity to benzo[d]isoxazole and other well characterized series and FEP calculations were performed for several tautomers of the deprotonated and protonated compounds since all these forms are potentially present owing to the experimental pKa values of representative compounds in the series. Deprotonated compounds are proposed to be the most important bound species owing to the significantly better agreement between their calculated and measured affinities compared to the protonated forms. FEP calculations were also used for the prediction of the affinities of compounds not previously tested as DAAO inhibitors and for a comparative structure–activity relationship study of the benzo[d]isoxazole and indazole series. Selected indazole derivatives were synthesized and their measured binding affinity towards DAAO was in good agreement with FEP predictions.

Estimating the correct binding modes of ligands in protein-ligand complexes is not only crucial in the drug discovery process, but also for elucidating potential toxicity mechanisms. In the current paper, we discuss and demonstrate a computational modelling protocol using the combination of docking, classical (cMD) and accelerated (aMD) molecular dynamics and free energy perturbation (FEP+ protocol) for identification of the binding modes of selected perfluorocarboxyl acids (PFCAs) in the PPARγ nuclear receptor. Initially, we employed both the regular and induced fit docking which failed to correctly predict the ligand binding modes and rank the compounds with respect to experimental free energies of binding, when they were docked into non-native X-ray structure. The cMD and aMD simulations identified the presence of multiple binding modes for these compounds, and the shorter chain PFCAs (C6-C8) continuously moved between a few energetically favourable binding conformations. These results demonstrate that the docking scoring function cannot rank compounds precisely in such cases, not due to its insufficiency, but because of the use of incorrect or only one unique bindings pose, neglecting the protein dynamics. Finally, based on MD predictions of binding conformations, the FEP+ sampling protocol was extended and then accurately reproduced experimental differences in the free energies. Thus, the preliminary MD simulations can also provide helpful information about correct set-up of the FEP+ calculations. These results show that the PFCAs binding modes were accurately predicted and the FEP+ protocol can be used to estimate free energies of binding of flexible molecules outside of typical drug-like compounds. Our in silico workflow revealed the main characteristics of the PFCAs, which are week PPARγ partial agonists and illustrated the importance of specific ligand-residue interactions within the LBD. This work also suggests a common workflow for identification of ligand binding modes, ligand-protein dynamics description and relative free energy calculations.

Transition state search is at the center of multiple types of computational chemical predictions related to mechanistic investigations, reactivity and regioselectivity predictions, and catalyst design. The process of finding transition states in practice is, however, a laborious multistep operation that requires significant user involvement. Here, we report a highly automated workflow designed to locate transition states for a given elementary reaction with minimal setup overhead. The only essential inputs required from the user are the structures of the separated reactants and products. The seamless workflow combining computational technologies from the fields of cheminformatics, molecular mechanics, and quantum chemistry automatically finds the most probable correspondence between the atoms in the reactants and the products, generates a transition state guess, launches a transition state search through a combined approach involving the relaxing string method and the quadratic synchronous transit, and finally validates the transition state via the analysis of the reactive chemical bonds and imaginary vibrational frequencies as well as by the intrinsic reaction coordinate method. Our approach does not target any specific reaction type, nor does it depend on training data; instead, it is meant to be of general applicability for a wide variety of reaction types. The workflow is highly flexible, permitting modifications such as a choice of accuracy, level of theory, basis set, or solvation treatment. Successfully located transition states can be used for setting up transition state guesses in related reactions, saving computational time and increasing the probability of success. The utility and performance of the method are demonstrated in applications to transition state searches in reactions typical for organic chemistry, medicinal chemistry, and homogeneous catalysis research. In particular, applications of our code to Michael additions, hydrogen abstractions, Diels-Alder cycloadditions, carbene insertions, and an enzyme reaction model involving a molybdenum complex are shown and discussed.

The emergence of multidrug-resistant Mycobacterium tuberculosis (Mtb) strains highlights the need to develop more efficacious and potent drugs. However, this goal is dependent on a comprehensive understanding of Mtb virulence protein effectors at the molecular level. Here, we used a post-expression cysteine (Cys)-to-dehydrolanine (Dha) chemical editing strategy to identify a water-mediated motif that modulates accessibility of the protein tyrosine phosphatase A (PtpA) catalytic pocket. Importantly, this water-mediated Cys-Cys non-covalent motif is also present in the phosphatase SptpA from Staphylococcus aureus, which suggests a potentially preserved structural feature among bacterial tyrosine phosphatases. The identification of this structural water provides insight into the known resistance of Mtb PtpA to the oxidative conditions that prevail within an infected host macrophage. This strategy could be applied to extend the understanding of the dynamics and function(s) of proteins in their native state and ultimately aid in the design of small-molecule modulators.

A series of acylguanidine beta secretase 1 (BACE1) inhibitors with modified scaffold and P3 pocket substituent was synthesized and studied with free energy perturbation (FEP) calculations. The resulting molecules showed potencies in enzymatic BACE1 inhibition assays up to 1 nM. The correlation between the predicted activity from the FEP calculations and the experimental activity was good for the P3 pocket substituents. The average mean unsigned error (MUE) between prediction and experiment was 0.68 ± 0.17 kcal/mol for the default 5 ns lambda window simulation time improving to 0.35 ± 0.13 kcal/mol for 40 ns. FEP calculations for the P2' pocket substituents on the same acylguanidine scaffold also showed good agreement with experiment and the results remained stable with repeated simulations and increased simulation time. It proved more difficult to use FEP calculations to study the scaffold modification from increasing 5 to 6 and 7 membered-rings. Although prediction and experiment were in agreement for short 2 ns simulations, as the simulation time increased the results diverged. This was improved by the use of a newly developed "Core Hopping FEP+" approach, which also showed improved stability in repeat calculations. The origins of these differences along with the value of repeat and longer simulation times are discussed. This work provides a further example of the use of FEP as a computational tool for molecular design.

Protein side-chain mutation is fundamental both to natural evolutionary processes and to the engineering of protein therapeutics, which constitute an increasing fraction of important medications. Molecular simulation enables the prediction of the effects of mutation on properties such as binding affinity, secondary and tertiary structure, conformational dynamics, and thermal stability. A number of widely differing approaches have been applied to these predictions, including sequence-based algorithms, knowledge-based potential functions, and all-atom molecular mechanics calculations. Free energy perturbation theory, employing all-atom and explicit-solvent molecular dynamics simulations, is a rigorous physics-based approach for calculating thermodynamic effects of, for example, protein side-chain mutations. Over the past several years, we have initiated an investigation of the ability of our most recent free energy perturbation methodology to model the thermodynamics of protein mutation for two specific problems: protein–protein binding affinities and protein thermal stability. We highlight recent advances in the field and outline current and future challenges.

This presentation will describe how structural biology, molecular pharmacology, and medicinal chemistry studies can be combined with molecular modeling and chemoinformatics analyses for a more accurate description and prediction of structural determinants of protein-ligand binding, functional activity, and selectivity.The challenges and possibilities of structural chemogenomics studies will be discussed, including the integration of large volumes of heterogeneous pharmacological and chemical data for different protein targets and the development of structure-based virtual screening and computer-aided drug design approaches to discover novel small molecule ligands with well defined functional activity and protein selectivity profiles.The potential of molecular dynamics simulation methods to complement hybrid structural biology studies will be demonstrated for the investigation the mechanisms of conformational selection and protein-ligand binding kinetics.In the final part of the presentation structural protein-ligand interaction databases will be described that link structure-based protein-ligand interaction maps to protein ligand topology and can be used as structural chemogenomics tools to navigate medicinal chemistry space.

The stability of folded proteins is critical to their biological function and for the efficacy of protein therapeutics. Predicting the energetic effects of protein mutations can improve our fundamental understanding of structural biology, the molecular basis of diseases, and possible routes to addressing those diseases with biological drugs. Identifying the effect of single amino acid point mutations on the thermodynamic equilibrium between the folded and unfolded states of a protein can pinpoint residues of critical importance that should be avoided in the process of improving other properties (affinity, solubility, viscosity, etc.) and suggest changes at other positions for increasing stability in protein engineering. Multiple computational tools have been developed for in silico predictions of protein stability in recent years, ranging from sequence-based empirical approaches to rigorous physics-based free energy methods. In this work, we show that FEP+, which is a free energy perturbation method based on all-atom molecular dynamics simulations, can provide accurate thermal stability predictions for a wide range of biologically relevant systems. Significantly, the FEP+ approach, while originally developed for relative binding free energies of small molecules to proteins and not specifically fitted for protein stability calculations, performs well compared to other methods that were fitted specifically to predict protein stability. Here, we present the broadest validation of a rigorous free energy-based approach applied to protein stability reported to date: 700+ single-point mutations spanning 10 different protein targets. Across the entire data set, we correctly classify the mutations as stabilizing or destabilizing in 84% of the cases, and obtain statistically significant predictions as compared with experiment [average error of ~1.6kcal/mol and coefficient of determination (R2) of 0.40]. This study demonstrates, for the first time in a large-scale validation, that rigorous free energy calculations can be used to predict changes in protein stability from point mutations without parameterization or system-specific customization, although further improvements should be possible with additional sampling and a better representation of the unfolded state of the protein. Here, we describe the FEP+ method as applied to protein stability calculations, summarize the large-scale retrospective validation results, and discuss limitations of the method, along with future directions for further improvements.