Alternative approaches to potential of mean force calculations: Free energy perturbation versus thermodynamic integration. Case study of some representative nonpolar interactions

We investigated the convergence behavior of potential of mean force (PMF) calculations using free energy perturbation (FEP), thermodynamic integration (TI), and “slow growth” (SG) techniques. The critical comparison of these alternative approaches is illustrated by the study of three different systems: two tagged argon atoms in a periodic box of argon, two methane molecules, and two benzene molecules maintained in a “T-shaped” conformation, both dimers embedded in a periodic box of water. The complete PMF simulations were carried out considering several protocols, in which the number of intermediate “λ” states, together with the amount of sampling per individual state, were varied. In most cases, as much as 1 ns of molecular dynamics (MD) sampling was used to derive each free energy profile. For the different systems examined, we find that FEP and TI unquestionably constitute robust computational methods leading to results of comparable accuracy. We also show that proper convergence of the free energy calculations, and further quantitative interpretation of the PMFs, requires total simulation times much higher than has been hitherto estimated. In some circumstances, the free energy profiles derived from FEP calculations tend to be slightly poorer than those obtained with TI, as a probable consequence of the greater sensitivity of FEP to the window spacing δλ. In the context of TI, and to a lesser extent FEP, simulations, it appears preferable to employ a limited number of “λ” points of the integrand involving extensive sampling, rather than numerous points with fewer samplings. Finally, we note that, at least in the case of nonpolar interactions, PMFs of reasonable quality can be generated using SG, and at a substantially lower cost than with either FEP or TI. © 1996 by John Wiley & Sons, Inc.     

Free energy difference calculations by thermodynamic integration: Difficulties in obtaining a precise value

Free energy difference calculations have been performed by the “slow growth” method of thermodynamic integration of the AMBER 3.0 molecular dynamics program for the mutation of a conformationally restricted threonine dipeptide, N-acetyl threonyl-N-methylamide, to the corresponding alanyl dipeptide. By varying the total simulation length, it has been determined that precise free energy values are obtained only for simulations of greater than 100 ps total simulation time length. By varying the starting configurations for simulations of the same length, it has been determined that averaging the free energies obtained from shorter simulations may not give precise answers. Possible reasons for this behavior are discussed.     

Free energy calculations,enhanced by a gaussian ansatz,for the “chemical work” distribution

The evaluation of the free energy is essential in molecular simulation because it is intimately related with the existence of multiphase equilibrium. Recently, it was demonstrated that it is possible to evaluate the Helmholtz free energy using a single statistical ensemble along an entire isotherm by accounting for the “chemical work” of transforming each molecule, from an interacting one, to an ideal gas. In this work, we show that it is possible to perform such a free energy perturbation over a liquid vapor phase transition. Furthermore, we investigate the link between a general free energy perturbation scheme and the novel nonequilibrium theories of Crook's and Jarzinsky. We find that for finite systems away from the thermodynamic limit the second law of thermodynamics will always be an inequality for isothermal free energy perturbations, resulting always to a dissipated work that may tend to zero only in the thermodynamic limit. The work, the heat, and the entropy produced during a thermodynamic free energy perturbation can be viewed in the context of the Crooks and Jarzinsky formalism, revealing that for a given value of the ensemble average of the “irreversible” work, the minimum entropy production corresponded to a Gaussian distribution for the histogram of the work. We propose the evaluation of the free energy difference in any free energy perturbation based scheme on the average irreversible “chemical work” minus the dissipated work that can be calculated from the variance of the distribution of the logarithm of the work histogram, within the Gaussian approximation. As a consequence, using the Gaussian ansatz for the distribution of the “chemical work,” accurate estimates for the chemical potential and the free energy of the system can be performed using much shorter simulations and avoiding the necessity of sampling the computational costly tails of the “chemical work.” For a more general free energy perturbation scheme that the Gaussian ansatz may not be valid, the free energy calculation can be expressed in terms of the moment generating function of the “chemical work” distribution. © 2014 Wiley Periodicals, Inc.     

Free‐energy differences between states with different conformational ensembles

Multiple conformations separated by high‐energy barriers represent a challenging problem in free‐energy calculations due to the difficulties in achieving adequate sampling. We present an application of thermodynamic integration (TI) in conjunction with the local elevation umbrella sampling (LE/US) method to improve convergence in alchemical free‐energy calculations. TI‐LE/US was applied to the guanosine triphosphate (GTP) to 8‐Br‐GTP perturbation, molecules that present high‐energy barriers between the anti and syn states and that have inverted preferences for those states. The convergence and reliability of TI‐LE/US was assessed by comparing with previous results using the enhanced‐sampling one‐step perturbation (OSP) method. A linear interpolation of the end‐state biasing potentials was sufficient to dramatically improve sampling along the chosen reaction coordinate. Conformational free‐energy differences were also computed for the syn and anti states and compared to experimental and theoretical results. Additionally, a coupled OSP with LE/US was carried out, allowing the calculation of conformational and alchemical free energies of GTP and 8‐substituted GTP analogs. © 2013 Wiley Periodicals, Inc.     

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.     

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     

Calculation of relative free energy differences for the covalent hydration of organic compounds: A combined quantum mechanical and free energy perturbation study

The relative free energy difference (ΔΔG_hyd) for the reversible addition of water to two unsaturated molecules is accurately computed using a combination of ab initio quantum mechanical calculations and free energy perturbation methods. Initial attempts to calculate the absolute hydration free energy difference (ΔG_hyd) for formaldehyde and trichloroacetaldehyde gave values that differed substantially from experimental results even after inclusion of electron correlation energy contributions using third-order (MP3) and fourth-order (MP4) Møller-Plesset perturbation theory and QCISD(T) correlation methods at the 6-31G** basis set level. Inaccuracies in ΔG_hyd were attributed to errors in the calculation of both ΔG_gas and ΔΔG_sol. Gas phase quantum mechanical free energies (ΔG_gas) varied significantly (2–3 kcal/mol) depending on the level of theory. Errors in ΔΔG_sol were attributed to slow convergence of the calculations using the thermodynamic cycle perturbation (TCP) method with explicit solvent. These errors were minimized or canceled, however, when relative hydration free energy differences (ΔΔG_hyd) were calculated using a combination of ab initio quantum mechanical calculations and free energy perturbation methods. Calculated values for a variety of aldehydes and ketones were consistent with experimental data. © 1995 John Wiley & Sons, Inc.     

Calculation of relative solvation free energy differences by thermodynamic perturbation method: Dependence of free energy results on simulation length

Molecular dynamics (MD) simulations in conjunction with the thermodynamic cycle perturbation approach has been used to calculate relative solvation free energies for acetone to acetaldehyde, acetone to pyruvic acid, acetone to 1,1,1-trifluoroacetone, acetone to 1,1,1-trichloroacetone, acetone to 2,3-butanedione, acetone to cyclopropanone, and formaldehyde hydrate to formaldehyde. To evaluate the dependence of relative solvation free energy convergence on MD simulation length and starting configuration two studies were performed. In the first study, each simulation started from the same well-equilibrated configuration and the length was varied from 153 to 1530 ps. In the second study, the relative solvation free energy differences were calculated starting from three different configurations and using 510 ps of MD simulation for each mutation. These results clearly indicate that, even for molecules with limited conformational flexibility, a simulation length of 510 ps or greater is required to obtain satisfactory convergence and, for the mutations of large structural changes between reactant and product, such as cyclopropanone to acetone, require much longer simulation lengths to achieve satisfactory convergence. These results also show that performing one long simulation is better than averaging results from three shortest simulations of the same length using different starting conformations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1018–1027, 1999     

Free Energy Calculations with Non-Equilibrium Methods: Applications of the Jarzynski Relationship

We present an introduction to the Jarzynski relationship that makes a strong connection, for a thermodynamic transformation, between the distribution of non-equilibrium work values and the corresponding equilibrium free energy differences. The relationship is discussed in the context of sampling issues, high level parallel computing and convergence criteria. We discuss three different applications by our group: mechanical unfolding of peptides, mixed quantum/classical free energy calculations in enzymes, and ligand escape pathways     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Valence energy in variational Monte Carlo: CuH dissociation energy

We show how to estimate the dissociation energy of CuH using the variational Monte Carlo method. The techniques involved are (i) an all-electron approach, (ii) a diffusion-only Metroplis algorithm which is well-suited for sampling the nodal regions properly, and (iii) a core-valence partitioning scheme such that the dissociation energy is estimated from the valence energies of CuH and Cu only. This approach avoids several of the approximations inherent in pseudopotential methods. Using relatively crude wave functions, we obtain an estimate of the dissociation energy and dipole moment with an accuracy on par with much more elaborate calculations in the literature. © 1996 John Wiley & Sons, Inc.     

Bennett's acceptance ratio and histogram analysis methods enhanced by umbrella sampling along a reaction coordinate in configurational space

Free energy perturbation, a method for computing the free energy difference between two states, is often combined with non-Boltzmann biased sampling techniques in order to accelerate the convergence of free energy calculations. Here we present a new extension of the Bennett acceptance ratio (BAR) method by combining it with umbrella sampling (US) along a reaction coordinate in configurational space. In this approach, which we call Bennett acceptance ratio with umbrella sampling (BAR-US), the conditional histogram of energy difference (a mapping of the 3N-dimensional configurational space via a reaction coordinate onto 1D energy difference space) is weighted for marginalization with the associated population density along a reaction coordinate computed by US. This procedure produces marginal histograms of energy difference, from forward and backward simulations, with higher overlap in energy difference space, rendering free energy difference estimations using BAR statistically more reliable. In addition to BAR-US, two histogram analysis methods, termed Bennett overlapping histograms with US (BOH-US) and Bennett-Hummer (linear) least square with US (BHLS-US), are employed as consistency and convergence checks for free energy difference estimation by BAR-US. The proposed methods (BAR-US, BOH-US, and BHLS-US) are applied to a 1-dimensional asymmetric model potential, as has been used previously to test free energy calculations from non-equilibrium processes. We then consider the more stringent test of a 1-dimensional strongly (but linearly) shifted harmonic oscillator, which exhibits no overlap between two states when sampled using unbiased Brownian dynamics. We find that the efficiency of the proposed methods is enhanced over the original Bennett's methods (BAR, BOH, and BHLS) through fast uniform sampling of energy difference space via US in configurational space. We apply the proposed methods to the calculation of the electrostatic contribution to the absolute solvation free energy (excess chemical potential) of water. We then address the controversial issue of ion selectivity in the K(+) ion channel, KcsA. We have calculated the relative binding affinity of K(+) over Na(+) within a binding site of the KcsA channel for which different, though adjacent, K(+) and Na(+) configurations exist, ideally suited to these US-enhanced methods. Our studies demonstrate that the significant improvements in free energy calculations obtained using the proposed methods can have serious consequences for elucidating biological mechanisms and for the interpretation of experimental data.     

Synergistic approach to improve "alchemical" free energy calculation in rugged energy surface

The authors present an integrated approach to "alchemical" free energy simulation, which permits efficient calculation of the free energy difference on rugged energy surface. The method is designed to obtain efficient canonical sampling for rapid free energy convergence. The proposal is motivated by the insight that both the exchange efficiency in the presently designed dual-topology alchemical Hamiltonian replica exchange method (HREM), and the confidence of the free energy determination using the overlap histogramming method, depend on the same criterion, viz., the overlaps of the energy difference histograms between all pairs of neighboring states. Hence, integrating these two techniques can produce a joint solution to the problems of the free energy convergence and conformational sampling in the free energy simulations, in which lambda parameter plays two roles to simultaneously facilitate the conformational sampling and improve the phase space overlap for the free energy determination. Specifically, in contrast with other alchemical HREM based free energy simulation methods, the dual-topology approach can ensure robust conformational sampling. Due to these features (a synergistic solution to the free energy convergence and canonical sampling, and the improvement of the sampling efficiency with the dual-topology treatment), the present approach, as demonstrated in the model studies of the authors, is highly efficient in obtaining accurate free energy differences, especially for the systems with rough energy landscapes.     

Accurate estimation of solvation free energy using polynomial fitting techniques

This report details an approach to improve the accuracy of free energy difference estimates using thermodynamic integration data (slope of the free energy with respect to the switching variable λ) and its application to calculating solvation free energy. The central idea is to utilize polynomial fitting schemes to approximate the thermodynamic integration data to improve the accuracy of the free energy difference estimates. Previously, we introduced the use of polynomial regression technique to fit thermodynamic integration data (Shyu and Ytreberg, J Comput Chem, 2009, 30, 2297). In this report we introduce polynomial and spline interpolation techniques. Two systems with analytically solvable relative free energies are used to test the accuracy of the interpolation approach. We also use both interpolation and regression methods to determine a small molecule solvation free energy. Our simulations show that, using such polynomial techniques and nonequidistant λ values, the solvation free energy can be estimated with high accuracy without using soft‐core scaling and separate simulations for Lennard‐Jones and partial charges. The results from our study suggest that these polynomial techniques, especially with use of nonequidistant λ values, improve the accuracy for ΔF estimates without demanding additional simulations. We also provide general guidelines for use of polynomial fitting to estimate free energy. To allow researchers to immediately utilize these methods, free software and documentation is provided via http://www.phys.uidaho.edu/ytreberg/software . © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Short-range potential functions in computer simulations of water and aqueous solutions

A review and comparative analyses of methods for restricting the range of molecular interactions within the concept of atom-atom potentials are presented. Emphasis is placed on the problem of calculating the electrostatic energy in models with periodic boundary conditions. Numerous calculations of the thermodynamic and structural characteristics of water using parallel Monte Carlo computations have shown that the use of functional forms simulating the electric potentials of “screened charges” provides very good results.     

Some practical aspects of free energy calculations from molecular dynamics simulation

In this work, we address two critical aspects of calculation of the free energy differences in molecular systems from molecular simulations. The first aspect involves checking whether the calculated free energy difference depends significantly on the extent of perturbation used for accomplishment of a given transformation. The second aspect of interest is to verify if the sampling errors in calculating the free energy differences between the wild-type molecule and a mutated one in its free state and in a complex are similar, or not, for a finite-length dynamic simulation. The reliability of the free energy estimates obtained from molecular simulations using thermodynamic cycles depends in part on this fact. For investigating these aspects, we use a self-transformation scheme in which a transformation of a part of a molecular system into itself is considered. We perform MD simulations of DNA fragments in which a part of a specific base is subjected to such a self-transformation. Results indicate that the estimated free energy differences do not depend significantly on the extent of perturbation used to achieve the transformation. Interestingly, the variation in the cumulative free energy difference, ΔA, with the coupling parameter, λ, depends significantly on the extent of perturbation. We examine the physical basis of the observed nature of the variation of the accumulated free energy difference, ΔA, against the λ value in the case of a self-transformation. In a thermodynamic cycle, the sampling errors due to the finite-length simulation for the molecular system are found to be similar to each other for the two perturbations (free and in a complex) justifying the use of such approach in calculating ΔΔA in molecular complexes. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 877–885, 1999     

On the direct calculation of the free energy of quantization for molecular systems in the condensed phase

Using the path integral formalism or the Feynman-Hibbs approach, various expressions for the free energy of quantization for a molecular system in the condensed phase can be derived. These lead to alternative methods to directly compute quantization free energies from molecular dynamics computer simulations, which were investigated with an eye to their practical use. For a test system of liquid neon, two methods are shown to be most efficient for a direct evaluation of the excess free energy of quantization. One of them makes use of path integral simulations in combination with a single-step free energy perturbation approach and was previously reported in the literature. The other method employs a Feynman-Hibbs effective Hamiltonian together with the thermodynamic integration formalism. However, both methods are found to give less accurate results for the excess free energy of quantization than the estimate obtained from explicit path integral calculations on the excess free energy of the neon liquid in the classical and quantum mechanical limit. Suggestions are made to make both methods more accurate.     

One- and multidimensional conformational free energy simulations

A new thermodynamic integration approach to conformational free energy simulations is presented. The method is applicable both to one-dimensional cases (reaction coordinates) and multidimensional situations (free energy surfaces). Analysis of the properties of the thermodynamic integration algorithm is used to formulate methods of calculating multidimensional free energy gradients. The method is applied to calculate the free energy profile for rotation around the central C—C bond of n-butane in the gas and liquid phase and to generate maps of the 18-dimensional free energy gradient with respect to all nine ϕ and nine ψ dihedrals of the decaalanine and deca-α-methylalanine peptides in vacuum. For n-butane essentially no change in the gauche–trans equilibrium between the gas and liquid is predicted within the CHARMM explicit hydrogen model, with the thermodynamic integration, thermodynamic perturbation, and direct simulation methods yielding free energy profiles that are identical within errors. For the decapeptides the right-handed helical region of conformational space is investigated. For decaalanine a minimum on the free energy surface is found in the vicinity of (ϕ, ψ) = (-64.5°, -42.5°) in the α-helix region; no minimum exists for 3₁₀-helix-type conformers. For deca-α-methylalanine free energy minima corresponding to both the α-helix at ( - 55.5°, - 51.5°) and the 3₁₀-helix at ( - 54°, - 29°) are found; the α-helix state is favored by about 4 kcal/mol and the barrier for the concerted 3₁₀-helix → α-helix transition is about 3 kcal/mol. The α-methylation also considerably increases the rigidity of the α-helix with respect to deformations. The computational efficiency, ease of generalization to calculations of multidimensional gradients, and analytical capability due to component analysis of free energy differences make the method a novel, powerful tool to improve the basic understanding of conformational equilibria of flexible molecules in condensed phases. A related scheme for energy minimization in the presence of holonomic constraints is also presented, allowing generation of adiabatic energy surfaces in constrained systems. © 1996 by John Wiley & Sons, Inc.     

Identifying low variance pathways for free energy calculations of molecular transformations in solution phase

Improving the efficiency of free energy calculations is important for many biological and materials design applications, such as protein-ligand binding affinities in drug design, partitioning between immiscible liquids, and determining molecular association in soft materials. We show that for any pair potential, moderately accurate estimation of the radial distribution function for a solute molecule is sufficient to accurately estimate the statistical variance of a sampling along a free energy pathway. This allows inexpensive analytical identification of low statistical error free energy pathways. We employ a variety of methods to estimate the radial distribution function (RDF) and find that the computationally cheap two-body "dilute gas" limit performs as well or better than 3D-RISM theory and other approximations for identifying low variance free energy pathways. With a RDF estimate in hand, we can search for pairwise interaction potentials that produce low variance. We give an example of a search minimizing statistical variance of solvation free energy over the entire parameter space of a generalized "soft core" potential. The free energy pathway arising from this optimization procedure has lower curvature in the variance and reduces the total variance by at least 50% compared to the traditional soft core solvation pathway. We also demonstrate that this optimized pathway allows free energies to be estimated with fewer intermediate states due to its low curvature. This free energy variance optimization technique is generalizable to solvation in any homogeneous fluid and for any type of pairwise potential and can be performed in minutes to hours, depending on the method used to estimate g(r).