Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

Alchemical free energy calculations provide a means for the accurate determination of free energies from atomistic simulations and are increasingly used as a tool for computational studies of protein-ligand interactions. Much attention has been placed on efficient ways to deal with the "endpoint singularity" effect that can cause simulation instabilities when changing the number of atoms. In this study we compare the performance of linear and several nonlinear transformation methods, among them separation shifted "soft core" scaling, for a popular test system, the hydration free energy of an amino acid side chain. All the nonlinear methods yield similar results if extensive sampling is performed, but soft core scaling provides smooth lambda curves that are best suited for commonly used numerical integration schemes. Additionally, results from a more flexible solute, hexane, will also be discussed.

We report AMBER force field parameters for biological simulations involving phosphorylation of serine, threonine, or tyrosine. The initial parameters used RESP fitting for the atomic partial charges and standard values for all other parameters such as Lennard-Jones coefficients. These were refined with the aid of a thermodynamic cycle consisting of experimentally determined pKa values, solvation energies from molecular dynamics free energy simulations, and gas phase basicities from QM calculations. A polarization energy term was included to account for the charge density change between the gas phase and solution, and solvation free energies were determined using thermodynamic integration. Parameter adjustment is required to obtain consistent thermodynamic results with better balanced electrostatic interactions between water and the phosphate oxygens. To achieve this, we modified the phosphate oxygen radii. A thermodynamically consistent parameter set can be derived for monoanions and requires an increase of the van der Waals phosphate oxygen radii of approximately 0.09 Å. Larger, residue-specific radii appear to be needed for dianions. The revised parameters developed here should be of particular interest for environments where simulations of multiple protonation states may be of interest.

Accurate ranking of compounds with regards to their binding affinity to a protein using computational methods is of great interest to pharmaceutical research. Physics-based free energy calculations are regarded as the most rigorous way to estimate binding affinity. In recent years, many retrospective studies carried out both in academia and industry have demonstrated its potential. Here, we present the results of large-scale prospective application of the FEP+ method in active drug discovery projects in an industry setting at Merck KGaA, Darmstadt, Germany. We compare these prospective data to results obtained on a new diverse, public benchmark of eight pharmaceutically relevant targets. Our results offer insights into the challenges faced when using free energy calculations in real-life drug discovery projects and identify limitations that could be tackled by future method development. The new public data set we provide to the community can support further method development and comparative benchmarking of free energy calculations.

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.

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.1 kcal/mol (95% CI [0.9, 1.3] kcal/mol) and correctly classifying mutations 0.9 as resistant or susceptible with 88% accuracy (95% CI [82,93]%). 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.