Improved convergence in dual-topology free energy calculations through use of harmonic restraints

When calculating free energy differences between two molecular systems by means of molecular dynamics simulation, accessory potential functions can help eliminate uninteresting configurational entropy contributions, improve convergence, and facilitate reversibility. In this work, we demonstrate that the use of a harmonic potential function to restrain key portions of a molecular system in a free energy perturbation dual-topology molecular dynamics approach dramatically improves convergence and precision of the calculation. Limitations of this technique are illustrated, and its use in conjunction with a fixed bond-length constraint is developed. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1278–1283, 1998     

Improving the efficiency and reliability of free energy perturbation calculations using overlap sampling methods

A challenge in free energy calculation for complex molecular systems by computer simulation is to obtain a reliable estimate within feasible computational time. In this study, we suggest an answer to this challenge by exploring a simple method, overlap sampling (OS), for producing reliable free-energy results in an efficient way. The formalism of the OS method is based on ensuring sampling of important overlapping phase space during perturbation calculations. This technique samples both forward and reverse free energy perturbation (FEP) to improve the free-energy calculation. It considers the asymmetry of the FEP calculation and features an ability to optimize both the precision and the accuracy of the measurement without affecting the simulation process itself. The OS method is tested at two optimization levels: no optimization (simple OS), and full optimization (equivalent to Bennett's method), and compared to conventional FEP techniques, including the widely used direct FEP averaging method, on three alchemical mutation systems: (a) an anion transformation in water solution, (b) mutation between methanol and ethane, and (c) alchemical change of an adenosine molecule. It is consistently shown that the reliability of free-energy estimates can be greatly improved using the OS techniques at both optimization levels, while the performance of Bennett's method is particularly striking. In addition, the efficiency of a calculation can be significantly improved because the method is able to (a) converge to the right answer quickly, and (b) work for large perturbations. The basic two-stage OS method can be extended to admit additional stages, if needed. We suggest that the OS method can be used as a general perturbation technique for computing free energy differences in molecular simulations.     

Free energy perturbation and molecular dynamics calculations of copper binding to azurin

Free energy perturbation/molecular dynamics simulations have been carried out on copper/azurin systems calculating the binding affinities of copper (II) ion to azurin either in the native or in the unfolded state. In order to test the validity of the strategy adopted for the calculations and to establish what force field is suitable for these kinds of calculations, three different force fields, AMBER, CVFF, and CFF, have been alternatively used for the calculations and the results have been compared with experimental data obtained by spectroscopic titrations of copper (II)/azurin solutions and denaturation experiments. Our findings have pointed out that only CFF gives satisfactory results, thus providing a reliable tool for copper binding simulations in copper protein.     

Basic ingredients of free energy calculations: A review

Methods to compute free energy differences between different states of a molecular system are reviewed with the aim of identifying their basic ingredients and their utility when applied in practice to biomolecular systems. A free energy calculation is comprised of three basic components: (i) a suitable model or Hamiltonian, (ii) a sampling protocol with which one can generate a representative ensemble of molecular configurations, and (iii) an estimator of the free energy difference itself. Alternative sampling protocols can be distinguished according to whether one or more states are to be sampled. In cases where only a single state is considered, six alternative techniques could be distinguished: (i) changing the dynamics, (ii) deforming the energy surface, (iii) extending the dimensionality, (iv) perturbing the forces, (v) reducing the number of degrees of freedom, and (vi) multi‐copy approaches. In cases where multiple states are to be sampled, the three primary techniques are staging, importance sampling, and adiabatic decoupling. Estimators of the free energy can be classified as global methods that either count the number of times a given state is sampled or use energy differences. Or, they can be classified as local methods that either make use of the force or are based on transition probabilities. Finally, this overview of the available techniques and how they can be best used in a practical context is aimed at helping the reader choose the most appropriate combination of approaches for the biomolecular system, Hamiltonian and free energy difference of interest. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Molecular dynamics simulations and free energy calculations on the enzyme 4‐hydroxyphenylpyruvate dioxygenase

4‐Hydroxyphenylpyruvate dioxygenase is a relevant target in both pharmaceutical and agricultural research. We report on molecular dynamics simulations and free energy calculations on this enzyme, in complex with 12 inhibitors for which experimental affinities were determined. We applied the thermodynamic integration approach and the more efficient one‐step perturbation. Even though simulations seem well converged and both methods show excellent agreement between them, the correlation with the experimental values remains poor. We investigate the effect of slight modifications on the charge distribution of these highly conjugated systems and find that accurate models can be obtained when using improved force field parameters. This study gives insight into the applicability of free energy methods and current limitations in force field parameterization. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

New approach to free energy of solvation applying continuum models to molecular dynamics simulation

A new approach to the calculation of the free energy of solvation from trajectories obtained by molecular dynamics simulation is presented. The free energy of solvation is computed as the sum of three contributions originated at the cavitation of the solute by the solvent, the solute-solvent nonpolar (repulsion and dispersion) interactions, and the electrostatic solvation of the solute. The electrostatic term is calculated based on ideas developed for the broadly used continuum models, the cavitational contribution from the excluded volume by the Claverie-Pierotti model, and the Van der Waals term directly from the molecular dynamics simulation. The proposed model is tested for diluted aqueous solutions of simple molecules containing a variety of chemically important functions: methanol, methylamine, water, methanethiol, and dichloromethane. These solutions were treated by molecular dynamics simulations using SPC/E water and the OPLS force field for the organic molecules. Obtained free energies of solvation are in very good agreement with experimental data.     

Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD‐based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907‐atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.     

Optimal designs for pairwise calculation: An application to free energy perturbation in minimizing prediction variability

Pairwise-based methods such as the free energy perturbation (FEP) method have been widely deployed to compute the binding free energy differences between two similar host–guest complexes. The calculated pairwise free energy difference is either directly adopted or transformed to absolute binding free energy for molecule rank ordering. We investigated, through both analytic derivations and simulations, how the selection of pairs in the experiment could impact the overall prediction precision. Our studies showed that (1) the estimated absolute binding free energy () derived from calculated pairwise differences (ΔΔG) through weighted least squares fitting is more precise in prediction than the pairwise difference values when the number of pairs is more than the number of ligands and (2) prediction precision is influenced by both the total number of pairs and the specifically selected pairs, the latter being critically important when the number of calculated pairs is limited. Furthermore, we applied optimal experimental design in pair selection and found that the optimally selected pairs can outperform randomly selected pairs in prediction precision. In an illustrative example, we showed that, upon weighing ligand structure similarity into design optimization, the weighted optimal designs are more efficient than the literature reported designs. This work provides a new approach to assess retrospective pairwise-based prediction results, and a method to design new prospective pairwise-based experiments for molecular lead optimization. © 2019 Wiley Periodicals, Inc.     

Modeling of novel CDK7 inhibitors activity by molecular dynamics and free energy perturbation methods

1. Download : Download high-res image (136KB)
2. Download : Download full-size image

     

Calculating zeros: Non-equilibrium free energy calculations

Free energy calculations on three model processes with theoretically known free energy changes have been performed using short simulation times. A comparison between equilibrium (thermodynamic integration) and non-equilibrium (fast growth) methods has been made in order to assess the accuracy and precision of these methods. The three processes have been chosen to represent processes often observed in biomolecular free energy calculations. They involve a redistribution of charges, the creation and annihilation of neutral particles and conformational changes. At very short overall simulation times, the thermodynamic integration approach using discrete steps is most accurate. More importantly, reasonable accuracy can be obtained using this method which seems independent of the overall simulation time. In cases where slow conformational changes play a role, fast growth simulations might have an advantage over discrete thermodynamic integration where sufficient sampling needs to be obtained at every λ-point, but only if the initial conformations do properly represent an equilibrium ensemble. From these three test cases practical lessons can be learned that will be applicable to biomolecular free energy calculations.     

Multipole electrostatics in hydration free energy calculations

Hydration free energy (HFE) is generally used for evaluating molecular solubility, which is an important property for pharmaceutical and chemical engineering processes. Accurately predicting HFE is also recognized as one fundamental capability of molecular mechanics force field. Here, we present a systematic investigation on HFE calculations with AMOEBA polarizable force field at various parameterization and simulation conditions. The HFEs of seven small organic molecules have been obtained alchemically using the Bennett Acceptance Ratio method. We have compared two approaches to derive the atomic multipoles from quantum mechanical calculations: one directly from the new distributed multipole analysis and the other involving fitting to the electrostatic potential around the molecules. Wave functions solved at the MP2 level with four basis sets (6-311G*, 6-311++G(2d,2p), cc-pVTZ, and aug-cc-pVTZ) are used to derive the atomic multipoles. HFEs from all four basis sets show a reasonable agreement with experimental data (root mean square error 0.63 kcal/mol for aug-cc-pVTZ). We conclude that aug-cc-pVTZ gives the best performance when used with AMOEBA, and 6-311++G(2d,2p) is comparable but more efficient for larger systems. The results suggest that the inclusion of diffuse basis functions is important for capturing intermolecular interactions. The effect of long-range correction to van der Waals interaction on the hydration free energies is about 0.1 kcal/mol when the cutoff is 12?, and increases linearly with the number of atoms in the solute/ligand. In addition, we also discussed the results from a hybrid approach that combines polarizable solute with fixed-charge water in the HFE calculation.     

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.     

Net charge changes in the calculation of relative ligand‐binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio‐)chemical thermodynamics. Many important endogenous receptor‐binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice‐summation scheme or a cutoff‐truncation scheme with Barker–Watts reaction‐field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest‐host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free‐energy calculation studies if changes in the net charge are involved. © 2013 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Efficient free energy calculations on small molecule host-guest systems - a combined linear interaction energy/one-step perturbation approach

Two efficient methods to calculate binding affinities of ligands with proteins have been critically evaluated by using sixteen small ligand host-guest complexes. It is shown that both the one-step (OS) perturbation method and the linear interaction energy (LIE) method have complementing strengths and weaknesses and can be optimally combined in a new manner. The OS method has a sound theoretical basis to address the free energy of cavity formation, whereas the LIE approach is more versatile and efficient to calculate the free energy of adding charges to such cavities. The off-term, which is neglected in the original LIE equation, can be calculated without additional costs from the OS, offering a powerful synergy between the two methods. The LIE/OS approach presented here combines the best of two worlds and for the model systems studied here, is more accurate than and as efficient as the original methods. It has a sound theoretical background and no longer requires any empirical parameters. The method appears very well suited for application in lead-optimization programmes in drug research, where the structure and dynamics of a series of molecules is of interest, together with an accurate calculation of the binding free energy.     

Partition coefficients of methylated DNA bases obtained from free energy calculations with molecular electron density derived atomic charges

Partition coefficients serve in various areas as pharmacology and environmental sciences to predict the hydrophobicity of different substances. Recently, they have also been used to address the accuracy of force fields for various organic compounds and specifically the methylated DNA bases. In this study, atomic charges were derived by different partitioning methods (Hirshfeld and Minimal Basis Iterative Stockholder) directly from the electron density obtained by electronic structure calculations in a vacuum, with an implicit solvation model or with explicit solvation taking the dynamics of the solute and the solvent into account. To test the ability of these charges to describe electrostatic interactions in force fields for condensed phases, the original atomic charges of the AMBER99 force field were replaced with the new atomic charges and combined with different solvent models to obtain the hydration and chloroform solvation free energies by molecular dynamics simulations. Chloroform–water partition coefficients derived from the obtained free energies were compared to experimental and previously reported values obtained with the GAFF or the AMBER‐99 force field. The results show that good agreement with experimental data is obtained when the polarization of the electron density by the solvent has been taken into account, and when the energy needed to polarize the electron density of the solute has been considered in the transfer free energy. These results were further confirmed by hydration free energies of polar and aromatic amino acid side chain analogs. Comparison of the two partitioning methods, Hirshfeld‐I and Minimal Basis Iterative Stockholder (MBIS), revealed some deficiencies in the Hirshfeld‐I method related to the unstable isolated anionic nitrogen pro‐atom used in the method. Hydration free energies and partitioning coefficients obtained with atomic charges from the MBIS partitioning method accounting for polarization by the implicit solvation model are in good agreement with the experimental values. © 2018 Wiley Periodicals, Inc.     

Improved precision and efficiency of free energy calculations for small systems using lambda-scaled atomic masses and separating conformational and transformational sampling

We present results showing the importance of appropriate treatment of atomic masses in molecular dynamics (MD)-based single topology free-energy perturbations (FEPs) on small molecule systems. The reversibility of gas phase simulations is significantly improved by scaling the atomic mass of mutated atoms with the lambda variable normally used for the scaling of energy terms. Because this effect is less pronounced for solvated systems, it will not cancel in estimates of the relative hydration free energy difference. The advantage of mass scaling is demonstrated by a null mutation of ethane to ethane and the calculation of the relative hydration free energy difference between ethane and n-propane. Furthermore, it is found that the simulation time necessary for converged MD/FEPs is prohibitively large for relative hydration free energy calculations on cyclic alkanes. Therefore, we explore an alternative free energy pathway including strongly constrained conformations to improve convergence in FEP simulations of flexible molecules.     

Systematic and statistical error in histogram-based free energy calculations

A common technique for the numerical calculation of free energies involves estimation of the probability density along a given coordinate from a set of configurations generated via simulation. The process requires discretization of one or more reaction coordinates to generate a histogram from which the continuous probability density is inferred. We show that the finite size of the intervals used to construct the histogram leads to quantifiable systematic error. The width of these intervals also determines the statistical error in the free energy, and the choice of the appropriate interval is therefore driven by the need to balance the two sources of error. We present a method for the construction of the optimal histogram for a given system, and show that the use of this technique requires little additional computational expense. We demonstrate the efficacy of the technique for a model system, and discuss how the principles governing the choice of discretization interval could be used to improve extended sampling techniques.     

Minimum MD simulation length required to achieve reliable results in free energy perturbation calculations: Case study of relative binding free energies of fructose‐1,6‐bisphosphatase inhibitors

In an attempt to establish the criteria for the length of simulation to achieve the desired convergence of free energy calculations, two studies were carried out on chosen complexes of FBPase‐AMP mimics. Calculations were performed for varied length of simulations and for different starting configurations using both conventional‐ and QM/MM‐FEP methods. The results demonstrate that for small perturbations, 1248 ps simulation time could be regarded a reasonable yardstick to achieve convergence of the results. As the simulation time is extended, the errors associated with free energy calculations also gradually tapers off. Moreover, when starting the simulation from different initial configurations of the systems, the results are not changed significantly, when performed for 1248 ps. This study carried on FBPase‐AMP mimics corroborates well with our previous successful demonstration of requirement of simulation time for solvation studies, both by conventional and ab initio FEP. The establishment of aforementioned criteria of simulation length serves a useful benchmark in drug design efforts using FEP methodologies, to draw a meaningful and unequivocal conclusion. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011