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

Abstract

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).

Guidelines for the analysis of free energy calculations

“…Specifically, the accuracy of TI is not a direct function of overlap in energy distributions but instead is a function of the average curvature $\langle \frac{{\partial }^{2}U}{\partial {\mathrm{\lambda }}^2}\rangle$ [30]. On the other hand, perturbation-based techniques do not depend on smoothness of the integrand, but rather on overlap in the sampled energy distributions.…”

Approaches for Calculating Solvation Free Energies and Enthalpies Demonstrated with an Update of the FreeSolv Database

“…98,102,103 A common soft-core form for Lennard-Jones potential between two particles i and j is: $$U_{ij}^{LJ}(r_{ij},\lambda )={\lambda }^{n}4{\epsilon }_{ij}\left(\frac{1}{{false[\alpha {false(1-\lambda false)}^m+{false(r_{ij}/{\sigma }_{ij}false)}^{6}false]}^2}-\frac{1}{\alpha {(1-\lambda )}^m+{(r_{ij}/{\sigma }_{ij})}^6}\right)$$ where ε ij and σ ij are the Lennard-Jones well-depth and lengthscale parameters, respectively, and α is a positive constant which should typically be set to 0.5. 103,104 The exponents m and n are most efficient at n = 1 and m = 1, but other values have been used too. 100,103–105 Improvements have been achieved by new soft-core functions that ease the problem with additional minima within the formulation of the original soft-core potential, 106 and alternate potentials that construct near optimal paths for alchemical simulations.…”

“…126 Our choice of soft-core is the so called 1-1-6 ( m and n equal to 1 in equation 8) which leads to statistical uncertainties approximately of the same size as uncertainties from simulations using optimized path soft-cores. 104 Pressure was maintained at 1 atm by the Parrinello-Rahman barostat. 127 Enthalpy and entropy decomposition required 60 nanosecond Langevin dynamics simulations, with two femtosecond timesteps at 298.15 K and 1 atm in water and in vacuo for each molecule in the database.…”

Nearest-neighbor parameter for inosine-cytosine pairs through a combined experimental and computational approach

“…Parameters associated with atoms (van der Waals radii and charges) were gradually changed between two end states (pairs A and B) to compute the free-energy difference. The 1-1-48 soft-core potential function (Pham and Shirts, 2011) was used to remove numerical instability during the calculation. Both charges and the van der Waals radii were treated with the soft-core potential.…”