Precise Binding Free Energy Calculations for Multiple Molecules Using an Optimal Measurement Network of Pairwise Differences

Rigorous physics-based methods to calculate binding free energies of protein-ligand complexes have become a valued component of structure-based drug design. Relative and absolute binding free energy calculations have been deployed prospectively in support of solving diverse drug discovery challenges. Here we review recent applications of binding free energy calculations to fragment growing and linking, scaffold hopping, binding pose validation, virtual screening, covalent enzyme inhibition, and positional analogue scanning. Furthermore, we discuss the merits of using protein models and highlight recent efforts to replace costly binding free energy calculations with predictions from machine learning models trained on a limited number of free energy perturbation or thermodynamic integration calculations thereby allowing for extended chemical space exploration.

We investigate the role of order/disorder transitions in alchemical simulations of protein-ligand absolute binding free energies. We show, in the context of a potential of mean force description, that for a benchmarking system (the complex of the L99A mutant of T4 lysozyme with 3-iodotoluene) and for a more challenging system relevant for medicinal applications (the complex of the farnesoid X receptor with inhibitor 26 from a recent D3R challenge) that order/disorder transitions can significantly hamper Hamiltonian replica exchange sampling efficiency and slow down the rate of equilibration of binding free energy estimates. We further show that our analytical model of alchemical binding combined with the formalism developed by Straub et al. for the treatment of order/disorder transitions of molecular systems can be successfully employed to analyze the transitions and help design alchemical schedules and soft-core functions that avoid or reduce the adverse effects of rare binding/unbinding transitions. The results of this work pave the way for the application of these techniques to the alchemical estimation with explicit solvation of hydration free energies and absolute binding free energies of systems undergoing order/disorder transitions.

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.

We propose a free energy calculation method for receptor-ligand binding, which have multiple binding poses that avoids exhaustive enumeration of the poses. For systems with multiple binding poses, the standard procedure is to enumerate orientations of the binding poses, restrain the ligand to each orientation, and then, calculate the binding free energies for each binding pose. In this study, we modify a part of the thermodynamic cycle in order to sample a broader conformational space of the ligand in the binding site. This modification leads to more accurate free energy calculation without performing separate free energy simulations for each binding pose. We applied our modification to simple model host-guest systems as a test, which have only two binding poses, by using a single decoupling method (SDM) in implicit solvent. The results showed that the binding free energies obtained from our method without knowing the two binding poses were in good agreement with the benchmark results obtained by explicit enumeration of the binding poses. Our method is applicable to other alchemical binding free energy calculation methods such as the double decoupling method (DDM) in explicit solvent. We performed a calculation for a protein-ligand system with explicit solvent using our modified thermodynamic path. The results of the free energy simulation along our modified path were in good agreement with the results of conventional DDM, which requires a separate binding free energy calculation for each of the binding poses of the example of phenol binding to T4 lysozyme in explicit solvent. © 2019 Wiley Periodicals, Inc.

Alchemical relative binding free energy (ΔΔG) calculations have shown high accuracy in predicting ligand binding affinity and have been used as important tools in computer-aided drug discovery and design. However, there has been limited research on the application of ΔΔG methods to membrane proteins despite the fact that these proteins represent a significant proportion of drug targets, play crucial roles in biological processes, and are implicated in numerous diseases. In this study, to predict the binding affinity of ligands to G protein-coupled receptors (GPCRs), we employed two ΔΔG calculation methods: thermodynamic integration (TI) with AMBER and the alchemical transfer method (AToM) with OpenMM. We calculated ΔΔG values for 53 transformations involving four class A GPCRs and evaluated the performance of AMBER-TI and AToM-OpenMM. In addition, we conducted tests using different numbers of windows and varying simulation times to achieve reliable ΔΔG results and to optimize resource utilization. Overall, both AMBER-TI and AToM-OpenMM show good agreement with the experimental data. Our results validate the applicability of AMBER-TI and AToM-OpenMM for optimization of lead compounds targeting membrane proteins.

Calculation of binding free energies, which characterize the stability and specificity of molecular interactions (e.g., drug formulations), has been established as an important application of modern computational chemistry. Though several methodologies currently exist within the literature, there is always a need for improvements in speed and accuracy that allow for efficient and accurate chemical screening. In this paper we explore the application of an unrestrained advanced sampling method, adaptive biasing force (ABF), to the calculation of standard binding free energies for a set of molecular host-guest systems for which the binding free energy is well-known: the cucurbit[7]uril-hosted systems utilized in the HYDROPHOBE challenge [ Assaf et al. J. Phys. Chem. B 2017 , 121 , 11144 - 11162 ]. We demonstrate that the use of the ABF advanced sampling method yields values systematically closer to experiment, with a lower uncertainty in the calculation, than thermodynamic integration and restrained sampling techniques commonly used in the literature. We hypothesize that this is due to the algorithm's ability to drive multiple binding and dissociation events through applied bias forces. These results show the promise of unrestrained advanced sampling methods along well-defined reaction coordinates to the calculation of binding energies in host-guest systems.

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.

We present an extension of the alchemical transfer method (ATM) for the estimation of relative binding free energies of molecular complexes applicable to conventional, as well as scaffold-hopping, alchemical transformations. Named ATM-RBFE, the method is implemented in the free and open-source OpenMM molecular simulation package and aims to provide a simpler and more generally applicable route to the calculation of relative binding free energies than what is currently available. ATM-RBFE is based on sound statistical mechanics theory and a novel coordinate perturbation scheme designed to swap the positions of a pair of ligands such that one is transferred from the bulk solvent to the receptor binding site while the other moves simultaneously in the opposite direction. The calculation is conducted directly in a single solvent box with a system prepared with conventional setup tools, without splitting of electrostatic and nonelectrostatic transformations, and without pairwise soft-core potentials. ATM-RBFE is validated here against the absolute binding free energies of the SAMPL8 GDCC host-guest benchmark set and against protein-ligand benchmark sets that include complexes of the estrogen receptor ERα and those of the methyltransferase EZH2. In each case the method yields self-consistent and converged relative binding free energy estimates in agreement with absolute binding free energies and reference literature values, as well as experimental measurements.

β-amyloid cleaving enzyme 1 (BACE1) is regarded as an important target of drug design toward the treatment of Alzheimer's disease (AD). In this study, three separate molecular dynamics (MD) simulations and calculations of binding free energies were carried out to comparatively determine the identification mechanism of BACE1 for three inhibitors, 60W, 954 and 60X. The analyses of MD trajectories indicated that the presence of three inhibitors influences the structural stability, flexibility and internal dynamics of BACE1. Binding free energies calculated by using solvated interaction energy (SIE) and molecular mechanics generalized Born surface area (MM-GBSA) methods reveal that the hydrophobic interactions provide decisive forces for inhibitor-BACE1 binding. The calculations of residue-based free energy decomposition suggest that the sidechains of residues L91, D93, S96, V130, Q134, W137, F169 and I179 play key roles in inhibitor-BACE1 binding, which provides a direction for future drug design toward the treatment of AD.

A semi-empirical method for estimation of binding free energy, recently proposed by Aqvist and coworkers, has been effectively tested in several protein-ligand binding cases. We have applied this linear interaction energy method to predict the binding of some N-benzy- loxycarbonyl-L-phenyl alanyl-L-alanine ketones with bovine cathepsin B and computed the respective absolute binding constants from averages of molecular dynamics simulations. It is found that the computer simulation results agree well with available experimental data and make it possible to understand better the origin of tight binding and inhibitor specificity of cathepsin B.

A method is proposed for the estimation of absolute binding free energy of interaction between proteins and ligands. The linear interaction energy method is combined with atom-bond electronegativity equalization method at σπ level Force field (fused into molecular mechanics) and generalized Born continuum model calculation of electrostatic solvation for the estimation of the absolute free energy of binding. The parameters of this method are calibrated by using a training set of 24 HIV-1 protease-inhibitor complexes (PDB entry 1AAQ). A correlation coefficient of 0.93 was obtained with a root mean square deviation of 0.70 kcal mol(-1) . This approach is further tested on seven inhibitor and protease complexes, and it provides small root mean square deviation between the calculated binding free energy and experimental binding free energy without reparametrization. By comparing the radii of gyration and the hydrogen bond distances between ligand and protein of three training model molecules, the consistent comparison result of binding free energy is obtained. It proves that this method of calculating the binding free energy with appropriate structural analysis can be applied to quickly assess new inhibitors of HIV-1 proteases. To test whether the parameters of this method can apply to other drug targets, we have validated this method for the drug target cyclooxygenase-2.

Effective virtual screening relies on our ability to make accurate prediction of protein-ligand binding, which remains a great challenge. In this work, utilizing the molecular-mechanics Poisson-Boltzmann (or Generalized Born) surface area approach, we have evaluated the binding affinity of a set of 156 ligands to seven families of proteins, trypsin β, thrombin α, cyclin-dependent kinase (CDK), cAMP-dependent kinase (PKA), urokinase-type plasminogen activator, β-glucosidase A, and coagulation factor Xa. The effect of protein dielectric constant in the implicit-solvent model on the binding free energy calculation is shown to be important. The statistical correlations between the binding energy calculated from the implicit-solvent approach and experimental free energy are in the range of 0.56-0.79 across all the families. This performance is better than that of typical docking programs especially given that the latter is directly trained using known binding data whereas the molecular mechanics is based on general physical parameters. Estimation of entropic contribution remains the barrier to accurate free energy calculation. We show that the traditional rigid rotor harmonic oscillator approximation is unable to improve the binding free energy prediction. Inclusion of conformational restriction seems to be promising but requires further investigation. On the other hand, our preliminary study suggests that implicit-solvent based alchemical perturbation, which offers explicit sampling of configuration entropy, can be a viable approach to significantly improve the prediction of binding free energy. Overall, the molecular mechanics approach has the potential for medium to high-throughput computational drug discovery.

Our interest is relative binding free energy (RBFE) calculations based on molecular simulations. These are promising tools for lead optimization in drug discovery, computing changes in binding free energy due to modifications of a lead compound. However, in the "alchemical" framework for RBFE calculations, some types of mutations have the potential to introduce error into the computed binding free energies. Here we explore the magnitude of this error in several different model binding calculations. We find that some of the calculations involving ring breaking have significant errors, and these errors are especially large in bridged ring systems. Since the error is a function of ligand strain, which is unpredictable in advance, we believe that ring breaking should be avoided when possible.

Qualitatively accurate calculations of ligand binding free energies (BFEs) as applied to chiral selective receptors using the MINTA (mode integration) algorithm are presented. Extensive conformational searches are first performed using a mixed mode Monte Carlo algorithm (MCMM/LMCS) followed by computation of the BFEs using the MINTA approach. The Merck (MMFFs) force field is used throughout. Deficiencies in the default MMFFs partial charges and torsion parameters assigned to atoms directly involved in substrate binding are systematically removed by reparametrization of these parameters using the program PAROPT and on the basis of Jaguar quantum mechanics calculations. The new force field parameters lead to a systematic improvement in qualitative trends of BFE differences (ΔΔGL-D) for different receptors as compared with the experimental enantioselectivities (ee's). The results suggest that this method can be extended to a larger and more complex class of receptors such as proteins.

Alchemical free energy (AFE) calculations can predictbinding affinitychanges as a function of structural modifications and have becomepowerful tools for lead optimization and drug discovery. Central tothe setup and performance of AFE calculations is the manner of mappingalchemical transformations, known as the topology model. Single, dual,and hybrid topology models have been used with various AFE methodsin the field. In recent works, λ-dynamics (λD) free energycalculations, specifically, have preferred the use of a hybrid multipletopology (HMT) for sampling multiple ligand perturbations. In thiswork, we evaluate a new topology method called ligand overlay (LO)for use with λD-based calculations, including the recently introducedλ-dynamics with a bias-updated Gibbs sampling (LaDyBUGS) approach.LO is a full multiple topology model that allows entire ligands tobe sampled and restrained within a λ-dynamics framework. Relativebinding free energies were computed with HMT or LO topology modelswith LaDyBUGS for 45 ligands across five protein benchmark systems.An overall Pearson R correlation of 0.98 and meanunsigned error of 0.32 kcal/mol were observed, suggesting that LOis a viable alternative topology model for λD-based calculations.We discuss the merits of using an HMT or LO model for future ligandstudies with λD or LaDyBUGS calculations.