Efficient free energy calculation of water across lipid membranes

An efficient free energy (FE) calculation of a water molecule to go across lipid membranes is presented. Both overlapping distribution and cavity insertion Widom methods are complementarily used. The former is useful for a dense region where water molecules are found, i.e., from the interfacial to bulk water region, while the latter works well in the low density region, i.e., the hydrocarbon region. Since both methods evaluate the excess chemical potential of water, the obtained FE profile is free from the fitting problem usually arisen when a FE difference method is used. A diphytanyl phosphatidylcholine bilayer is used for our test calculations. An excellent and fast convergence of the chemical potential is obtained when each method is applied for the appropriate region. The estimated FE barrier using the Ewald method for the electrostatic interaction is approximately 7.2 kcal/mol, which is higher than that using the interaction cutoff of 20 A by about 0.9 kcal/mol.     

Linear-scaling soft-core scheme for alchemical free energy calculations

Alchemical free energy calculations involving the removal or insertion of atoms into condensed phase systems generally make use of soft-core scaling of nonbonded interactions, designed to circumvent numerical instabilities that arise from weakly interacting "hard" atoms in close proximity. Current methods model soft-core atoms by introducing a nonlinear dependence between the shape of the interaction potential and the strength of the interaction. In this article, we propose a soft-core method that avoids introducing such a nonlinear dependence, through the application of a smooth flattening of the potential energy only in a region that is energetically accessible under normal conditions. We discuss the benefits that this entails and explore a selection of applications, including enhanced methods for the estimation of free energy differences and for the automated optimization of the placement of intermediate states in multistage alchemical calculations.     

Predicting relative binding affinities of non-peptide HIV protease inhibitors with free energy perturbation calculations

The relative binding free energies in HIV protease of haloperidol thioketal (THK) and three of its derivatives were examined with free energy calculations. THK is a weak inhibitor (IC50 = 15 M) for which two cocrystal structures with HIV type 1 proteases have been solved [Rutenber, E. et al., J. Biol. Chem., 268 (1993) 15343]. A THK derivative with a phenyl group on C2 of the piperidine ring was expected to be a poor inhibitor based on experiments with haloperidol ketal and its 2- phenyl derivative (Caldera, P., personal communication). Our calculations predict that a 5-phenyl THK derivative, suggested based on examination of the crystal structure, will bind significantly better than THK. Although there are large error bars as estimated from hysteresis, the calculations predict that the 5-phenyl substituent is clearly favored over the 2-phenyl derivative as well as the parent compound. The unfavorable free energies of solvation of both phenyl THK derivatives relative to the parent compound contributed to their predicted binding free energies. In a third simulation, the change in binding free energy for 5-benzyl THK relative to THK was calculated. Although this derivative has a lower free energy in the protein, its decreased free energy of solvation increases the predicted G(bind) to the same range as that of the 2-phenyl derivative.     

Rapid approximate calculation of water binding free energies in the whole hydration domain of (bio)macromolecules

The evaluation of water binding free energies around solute molecules is important for the thermodynamic characterization of hydration or association processes. Here, a rapid approximate method to estimate water binding free energies around (bio)macromolecules from a single molecular dynamics simulation is presented. The basic idea is that endpoint free‐energy calculation methods are applied and the endpoint quantities are monitored on a three‐dimensional grid around the solute. Thus, a gridded map of water binding free energies around the solute is obtained, that is, from a single short simulation, a map of favorable and unfavorable water binding sites can be constructed. Among the employed free‐energy calculation methods, approaches involving endpoint information pertaining to actual thermodynamic integration calculations or endpoint information as exploited in the linear interaction energy method were examined. The accuracy of the approximate approaches was evaluated on the hydration of a cage‐like molecule representing either a nonpolar, polar, or charged water binding site and on α‐ and β‐cyclodextrin molecules. Among the tested approaches, the linear interaction energy method is considered the most viable approach. Applying the linear interaction energy method on the grid around the solute, a semi‐quantitative thermodynamic characterization of hydration around the whole solute is obtained. Disadvantages are the approximate nature of the method and a limited flexibility of the solute. © 2016 Wiley Periodicals, Inc.     

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.     

Walking freely in the energy and temperature space by the modified replica exchange molecular dynamics method

Replica Exchange Molecular Dynamics (REMD) method is a powerful sampling tool in molecular simulations. Recently, we made a modification to the standard REMD method. It places some inactive replicas at different temperatures as well as the active replicas. The method completely decouples the number of the active replicas and the number of the temperature levels. In this article, we make a further modification to our previous method. It uses the inactive replicas in a different way. The inactive replicas first sample in their own knowledge‐based energy databases and then participate in the replica exchange operations in the REMD simulation. In fact, this method is a hybrid between the standard REMD method and the simulated tempering method. Using different active replicas, one can freely control the calculation quantity and the convergence speed of the simulation. To illustrate the performance of the method, we apply it to some small models. The distribution functions of the replicas in the energy space and temperature space show that the modified REMD method in this work can let the replicas walk freely in both of the two spaces. With the same number of the active replicas, the free energy surface in the simulation converges faster than the standard REMD. © 2016 Wiley Periodicals, Inc.     

On the Relative Merits of Equilibrium and Non‐Equilibrium Simulations for the Estimation of Free‐Energy Differences

The possibility of estimating equilibrium free‐energy profiles from multiple non‐equilibrium simulations using the fluctuation–dissipation theory or the relation proposed by Jarzynski has attracted much attention. Although the Jarzynski estimator has poor convergence properties for simulations far from equilibrium, corrections have been derived for cases in which the work is Gaussian distributed. Here, we examine the utility of corrections proposed by Gore and collaborators using a simple dissipative system as a test case. The system consists of a single methane‐like particle in explicit water. The Jarzynski equality is used to estimate the change in free energy associated with pulling the methane particle a distance of 3.9 nm at rates ranging from ～0.1 to 100 m s^?1. It is shown that although the corrections proposed by Gore and collaborators have excellent numerical performance, the profiles still converge slowly. Even when the corrections are applied in an ideal case where the work distribution is necessarily Gaussian, performing simulations under quasi‐equilibrium conditions is still most efficient. Furthermore, it is shown that even for a single methane molecule in water, pulling rates as low as 1 m s^?1 can be problematic. The implications of this finding for studies in which small molecules or even large biomolecules are pulled through inhomogeneous environments at similar pulling rates are discussed.     

Combining the biased and unbiased sampling strategy into one convenient free energy calculation method

Constructing a free energy landscape for a large molecule is difficult. One has to use either a high temperature or a strong driving force to enhance the sampling on the free energy barriers. In this work, we propose a mixed method that combines these two kinds of acceleration strategies into one simulation. First, it applies an adaptive biasing potential to some replicas of the molecule. These replicas are particularly accelerated in a collective variable space. Second, it places some unbiased and exchangeable replicas at various temperature levels. These replicas generate unbiased sampling data in the canonical ensemble. To improve the sampling efficiency, biased replicas transfer their state variables to the unbiased replicas after equilibrium by Monte Carlo trial moves. In comparison to previous integrated methods, it is more convenient for users. It does not need an initial reference biasing potential to guide the sampling of the molecule. And it is also unnecessary to insert many replicas for the requirement of passing the free energy barriers. The free energy calculation is accomplished in a single stage. It samples the data as fast as a biased simulation and it processes the data as simple as an unbiased simulation. The method provides a minimalist approach to the construction of the free energy landscape. © 2019 Wiley Periodicals, Inc.     

On the accurate estimation of free energies using the jarzynski equality

The Jarzynski equality is one of the most widely celebrated and scrutinized nonequilibrium work theorems, relating free energy to the external work performed in nonequilibrium transitions. In practice, the required ensemble average of the Boltzmann weights of infinite nonequilibrium transitions is estimated as a finite sample average, resulting in the so-called Jarzynski estimator, . Alternatively, the second-order approximation of the Jarzynski equality, though seldom invoked, is exact for Gaussian distributions and gives rise to the Fluctuation-Dissipation estimator . Here we derive the parametric maximum-likelihood estimator (MLE) of the free energy considering unidirectional work distributions belonging to Gaussian or Gamma families, and compare this estimator to . We further consider bidirectional work distributions belonging to the same families, and compare the corresponding bidirectional to the Bennett acceptance ratio () estimator. We show that, for Gaussian unidirectional work distributions, is in fact the parametric MLE of the free energy, and as such, the most efficient estimator for this statistical family. We observe that and perform better than and , for unidirectional and bidirectional distributions, respectively. These results illustrate that the characterization of the underlying work distribution permits an optimal use of the Jarzynski equality. © 2018 Wiley Periodicals, Inc.     

Comparison of thermodynamic integration and Bennett acceptance ratio for calculating relative protein‐ligand binding free energies

The performances of Bennett's acceptance ratio method and thermodynamic integration (TI) for the calculation of free energy differences in protein simulations are compared. For the latter, the standard trapezoidal rule, Simpson's rule, and Clenshaw‐Curtis integration are used as numerical integration methods. We evaluate the influence of the number and definition of intermediate states on the precision, accuracy, and efficiency of the free energy calculations. Our results show that non‐equidistantly spaced intermediate states are in some cases beneficial for the TI methods. Using several combinations of softness parameters and the λ power dependence, it is shown that these benefits are strongly dependent on the shape of the integrand. Although TI is more user‐friendly due to its simplicity, it was found that Bennett's acceptance ratio method is the more efficient method. It is also the least dependent on the choice of the intermediate states, making it more robust than TI. © 2013 Wiley Periodicals, Inc.     

Constructing a multidimensional free energy surface like a spider weaving a web

Complete free energy surface in the collective variable space provides important information of the reaction mechanisms of the molecules. But, sufficient sampling in the collective variable space is not easy. The space expands quickly with the number of the collective variables. To solve the problem, many methods utilize artificial biasing potentials to flatten out the original free energy surface of the molecule in the simulation. Their performances are sensitive to the definitions of the biasing potentials. Fast‐growing biasing potential accelerates the sampling speed but decreases the accuracy of the free energy result. Slow‐growing biasing potential gives an optimized result but needs more simulation time. In this article, we propose an alternative method. It adds the biasing potential to a representative point of the molecule in the collective variable space to improve the conformational sampling. And the free energy surface is calculated from the free energy gradient in the constrained simulation, not given by the negative of the biasing potential as previous methods. So the presented method does not require the biasing potential to remove all the barriers and basins on the free energy surface exactly. Practical applications show that the method in this work is able to produce the accurate free energy surfaces for different molecules in a short time period. The free energy errors are small in the cases of various biasing potentials. © 2017 Wiley Periodicals, Inc.     

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     

Molecular dynamics simulation of water near nanostructured hydrophobic surfaces: interfacial energies.

We present results from molecular dynamics simulations of water near structured hydrophobic surfaces. The surface structures reported herein are a planar alkane crystal as a reference and crystals with a hole and a protrusion of approximately 2.5 nm diameter and 0.5 nm depth or height. All indicators show that surface structuring increases the hydrophobicity: The water density is reduced near the structure elements, and the number of residual contacts between water and the surface decreases by about 40 % with respect to the planar surface. Thermodynamic integration shows that the interfacial energy of the structured surfaces is about 7 mJ m(-2) higher for structured surfaces than for the planar surface. The hydrophobicity increases by a similar amount for the hole and the protrusion geometries compared to the planar surface.     

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.     

Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies

A free energy perturbation (FEP) method was developed that uses ab initio quantum mechanics (QM) for treating the solute molecules and molecular mechanics (MM) for treating the surroundings. Like our earlier results using AM1 semi empirical QMs, the ab initio QM/MM-based FEP method was shown to accurately calculate relative solvation free energies for a diverse set of small molecules that differ significantly in structure, aromaticity, hydrogen bonding potential, and electron density. Accuracy was similar to or better than conventional FEP methods. The QM/MM-based methods eliminate the need for time-consuming development of MM force field parameters, which are frequently required for drug-like molecules containing structural motifs not adequately described by MM. Future automation of the method and parallelization of the code for Linux 128/256/512 clusters is expected to enhance the speed and increase its use for drug design and lead optimization.     

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.