Overlap perturbation methods for computing alchemical free energy changes: variants, generalizations and evaluations

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.

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.

Guanine to inosine (G → I) substitution has often been used to study various properties of nucleic acids. Inosine differs from guanine only by loss of the N2 amino group, while both bases have similar electrostatic potentials. Therefore, G → I substitution appears to be optimally suited to probe structural and thermodynamics effects of single H-bonds and atomic groups. However, recent experiments have revealed substantial difference in free energy impact of G → I substitution in the context of B-DNA and A-RNA canonical helices, suggesting that the free energy changes reflect context-dependent balance of energy contributions rather than intrinsic strength of a single H-bond. In the present study, we complement the experiments by free energy computations using thermodynamics integration method based on extended explicit solvent molecular dynamics simulations. The computations successfully reproduce the basic qualitative difference in free energy impact of G → I substitution in B-DNA and A-RNA helices although the magnitude of the effect is somewhat underestimated. The computations, however, do not reproduce the salt dependence of the free energy changes. We tentatively suggest that the different effect of G → I substitution in A-RNA and B-DNA may be related to different topologies of these helices, which affect the electrostatic interactions between the base pairs and the negatively charged backbone. Limitations of the computations are briefly discussed.

Stressing crystals with solutes: Effects of added solutes on crystalline caffeine and their relevance to determining transfer free energies

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.

Many biomolecules undergo conformational changes associated with allostery or ligand binding. Observing these changes in computer simulations is difficult if their timescales are long. These calculations can be accelerated by observing the transition on an auxiliary free energy surface with a simpler Hamiltonian and connecting this free energy surface to the target free energy surface with free energy calculations. Here, we show that the free energy legs of the cycle can be replaced with energy representation (ER) density functional approximations. We compute: (1) The conformational free energy changes for alanine dipeptide transitioning from the right-handed free energy basin to the left-handed basin and (2) the free energy difference between the open and closed conformations of β-cyclodextrin, a "host" molecule that serves as a model for molecular recognition in host-guest binding. β-cyclodextrin contains 147 atoms compared to 22 atoms for alanine dipeptide, making β-cyclodextrin a large molecule for which to compute solvation free energies by free energy perturbation or integration methods and the largest system for which the ER method has been compared to exact free energy methods. The ER method replaced the 28 simulations to compute each coupling free energy with two endpoint simulations, reducing the computational time for the alanine dipeptide calculation by about 70% and for the β-cyclodextrin by > 95%. The method works even when the distribution of conformations on the auxiliary free energy surface differs substantially from that on the target free energy surface, although some degree of overlap between the two surfaces is required. © 2016 Wiley Periodicals, Inc.

The cytoplasm of a cell is crowded with many different kinds of macromolecules. The macromolecular crowding affects the thermodynamics and kinetics of biological reactions in a living cell, such as protein folding, association, and diffusion. Theoretical and simulation studies using simplified models focus on the essential features of the crowding effects and provide a basis for analyzing experimental data. In most of the previous studies on the crowding effects, a uniform crowder size is assumed, which is in contrast to the inhomogeneous size distribution of macromolecules in a living cell. Here, we evaluate the free energy changes upon macromolecular association in a cell-like inhomogeneous crowding system via a theory of hard-sphere fluids and free energy calculations using Brownian dynamics trajectories. The inhomogeneous crowding model based on 41 different types of macromolecules represented by spheres with different radii mimics the physiological concentrations of macromolecules in the cytoplasm of Mycoplasma genitalium. The free energy changes of macromolecular association evaluated by the theory and simulations were in good agreement with each other. The crowder size distribution affects both specific and nonspecific molecular associations, suggesting that not only the volume fraction but also the size distribution of macromolecules are important factors for evaluating in vivo crowding effects. This study relates in vitro experiments on macromolecular crowding to in vivo crowding effects by using the theory of hard-sphere fluids with crowder-size heterogeneity.

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.

A hybrid quantum and molecular mechanical (QM/MM) free energy perturbation (FEP) method is implemented in the context of molecular dynamics (MD). The semiempirical quantum mechanical (QM) Hamiltonian (Austin Model 1) represents solute molecules, and the molecular mechanical (MM) CHARMM force field describes the water solvent. The QM/MM FEP method is used to calculate the free energy changes in aqueous solution for (1) a proton transfer from methanol to imidazole and (2) a hydride transfer from methoxide to nicotinamide. The QM/MM interaction energies between the solute and solvent arc calibrated to emulate the solute-solvent interaction energies determined at the Hartee-Fock 6-31G(d) level of ab initio theory. The free energy changes for the proton and hydride transfers are calculated to be 15.1 and {minus}6.3 kcal/mol, respectively, which compare favorably with the corresponding experimental values of 12.9 and {minus}7.4 kcal/mol. An estimate of the reliability of the calculations is obtained through the computation of the forward (15.1 and {minus}6.3 kcal/mol) and backward ({minus}14.1 and 9.1 kcal/mol)free energy changes. The reasonable correspondence between these two independent calculations suggests that adequate phase space sampling is obtained along the reaction pathways chosen to transform the proton and hydride systems between their respective reactant and product states.

Rabbit anti-dinitrophenyl antibody from a serum pool was obtained as five fractions of purified specific antibody by a limiting antigen precipitation method. Each fraction had a different binding affinity for epsilon-N-2,4-dinitrophenyl-L-lysine. The free energy changes for hapten binding to the five antibody fractions varied from -8.35 to -10.0 kcal/mol. An average deltaH of -13.9 kcal/mol was measured for the fractions with a batch calorimeter. The results indicate no significant correlation between enthalpy changes and free energy changes. However, a statistically significant correlation between the free energy changes and the entropy changes was found. The enthalpy of the anti-dinitrophenyl antibody interaction with multivalent dinitrophenyl human serum albumin was determined. These are the first enthalpy measurements of an antibody antigen reaction in which the intrinsic binding enthalpy between the antibody and the determinant group is known. The deltaH for the antigen binding reaction was -10.1 kcal/mol which is 3.8 kcal/mol less exothermic than the deltaH for the hapten binding reaction. The interactions that could lead to such a difference in enthalpy are discussed.

A new practical approach to studying enzyme reactions by combining ab initio QM/MM calculations with free energy perturbation is presented. An efficient iterative optimization procedure has been developed to determine optimized structures and minimum energy paths for a system with thousands of atoms on the ab initio QM/MM potential: the small QM sub-system is optimized using a quasi-Newton minimizer in redundant internal coordinates with ab initio QM/MM calculations, while the large MM sub-system is minimized by the truncated Newton method in Cartesian coordinates with only molecular mechanical calculations. The above two optimization procedures are performed iteratively until they converge. With the determined minimum energy paths, free energy perturbation calculations are carried out to determine the change in free energy along the reaction coordinate. Critical to the success of the iterative optimization procedure and the free energy calculations is the smooth connection between the QM and MM regions provided by a recently proposed pseudobond QM/MM approach [J. Chem. Phys. 110, 46 (1999)]. The methods have been demonstrated by studying the initial proton transfer step in the reaction catalyzed by the enzyme triosephosphate isomerase (TIM).

We have ranked 13 different TIBO derivatives with respect to their relative free energies of binding using two approximate computational methods: adaptive chemical Monte Carlo/molecular dynamics (CMC/MD) and Poisson-Boltzmann/solvent accessibility (PB/SA) calculations. Eight of these derivatives have experimentally determined binding affinities. The remaining new derivatives were constructed based on contour maps around R86183 (8Cl-TIBO), generated with the program PROFEC (pictorial representation of free energy changes). The rank order among the derivatives with known binding affinity was in good agreement with experimental results for both methods, with average errors in the binding free energies of 1. 0 kcal/mol for CMC/MD and 1.3 kcal/mol for the PB/SA method. With both methods, we found that one of the new derivatives was predicted to bind 1-2 kcal/mol better than R86183, which is the hitherto most tightly binding derivative. This result was subsequently supported by the most rigorous free energy computational methods: free energy perturbation (FEP) and thermodynamic integration (TI). The strategy we have used here should be generally useful in structure-based drug optimization. An initial ligand is derivatized based on PROFEC suggestions, and the derivatives are ranked with CMC/MD and PB/SA to identify promising compounds. Since these two methods rely on different sets of approximations, they serve as a good complement to each other. Predictions of the improved affinity can be reinforced with FEP or TI and the best compounds synthesized and tested. Such a computational strategy would allow many different derivatives to be tested in a reasonable time, focusing synthetic efforts on the most promising modifications.

In the context of the SAMPL5 blinded challenge standard free energies of binding were predicted for a dataset of 22 small guest molecules and three different host molecules octa-acids (OAH and OAMe) and a cucurbituril (CBC). Three sets of predictions were submitted, each based on different variations of classical molecular dynamics alchemical free energy calculation protocols based on the double annihilation method. The first model (model A) yields a free energy of binding based on computed free energy changes in solvated and host-guest complex phases; the second (model B) adds long range dispersion corrections to the previous result; the third (model C) uses an additional standard state correction term to account for the use of distance restraints during the molecular dynamics simulations. Model C performs the best in terms of mean unsigned error for all guests (MUE [Formula: see text]-95% confidence interval) for the whole data set and in particular for the octa-acid systems (MUE [Formula: see text]). The overall correlation with experimental data for all models is encouraging ([Formula: see text]). The correlation between experimental and computational free energy of binding ranks as one of the highest with respect to other entries in the challenge. Nonetheless the large MUE for the best performing model highlights systematic errors, and submissions from other groups fared better with respect to this metric.

Molecular dynamics (MD) simulation combined with free energy perturbation (FEP) methods have been used to study the key structural differences and relative free energies for the binding of 6-methyl-N5-deazapterin (N8 protonated) and the 8-substituted compound, 6,8-dimethyl-N5-deazapterin (N3 protonated), to dihydrofolate reductase (DHFR). The free energy changes have been calculated using a variety of initial X-ray coordinates derived from bacterial and vertebrate (including human) DHFRs, and both with and without the reduced cofactor nicotinamide adenine dinucleotide (NADPH) bound. Given a sufficiently long simulation time for the FEP calculations (ca. 200 ps), all structures obtained after mutating 6,8-methyl-N5-deazapterin to 6-methyl-N5-deazapterin exhibited hydrogen bond formation between a backbone carbonyl group of DHFR and H(N8) of 6-methyl-N5-deazapterin, analogous to that found in the X-ray crystal structure of N5-deazafolate(N8 protonated) bound to human DHFR. However, both simulation and experiment suggest this additional H-bonding does not greatly enhance thermodynamic stability, with experiment indicating at most a factor of 2 difference in the relative affinities of the two ligand cations for vertebrate DHFR. Moreover, a binding differential of 10 in favour of the protonated 8-substituted compound is found experimentally for bacterial DHFR. The MD/FEP calculations suggest that the relative cost of ligand desolvation may largely cancel the lowering of free energy obtained in the active site, resulting in predicted binding differences within the range indicated by the vertebrate and bacterial DHFR experiments. However, the theoretical free energy changes could not be obtained with the accuracy required for the rationalization of the observed species dependence. While sampling difficulties are known to be inherent in MD simulation methodologies, these studies with several initial coordinate sets have demonstrated the contribution of coordinate choice to this problem. The results indicate that for demanding protein-ligand binding problems such as this one, the accuracy of the method may be no better than ± 2 kcal/mol.

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.