Overcoming free energy barriers using unconstrained molecular dynamics simulations

Association of unconstrained molecular dynamics (MD) and the formalisms of thermodynamic integration and average force [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)] have been employed to determine potentials of mean force. When implemented in a general MD code, the additional computational effort, compared to other standard, unconstrained simulations, is marginal. The force acting along a chosen reaction coordinate xi is estimated from the individual forces exerted on the chemical system and accumulated as the simulation progresses. The estimated free energy derivative computed for small intervals of xi is canceled by an adaptive bias to overcome the barriers of the free energy landscape. Evolution of the system along the reaction coordinate is, thus, limited by its sole self-diffusion properties. The illustrative examples of the reversible unfolding of deca-L-alanine, the association of acetate and guanidinium ions in water, the dimerization of methane in water, and its transfer across the water liquid-vapor interface are examined to probe the efficiency of the method.     

Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes

A hybrid quantum/classical path integral Monte Carlo (QC-PIMC) method for calculating the quantum free energy barrier for hydrogen transfer reactions in condensed phases is presented. In this approach, the classical potential of mean force along a collective reaction coordinate is calculated using umbrella sampling techniques in conjunction with molecular dynamics trajectories propagated according to a mapping potential. The quantum contribution is determined for each configuration along the classical trajectory with path integral Monte Carlo calculations in which the beads move according to an effective mapping potential. This type of path integral calculation does not utilize the centroid constraint and can lead to more efficient sampling of the relevant region of conformational space than free-particle path integral sampling. The QC-PIMC method is computationally practical for large systems because the path integral sampling for the quantum nuclei is performed separately from the classical molecular dynamics sampling of the entire system. The utility of the QC-PIMC method is illustrated by an application to hydride transfer in the enzyme dihydrofolate reductase. A comparison of this method to the quantized classical path and grid-based methods for this system is presented.     

Recovering the Crooks equation for dynamical systems in the isothermal-isobaric ensemble: a strategy based on the equations of motion

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)] relates the work done on a system during a nonequilibrium transformation to the free energy difference between the final and the initial state of the transformation. Recently, the authors have derived the Crooks equation for systems in the canonical ensemble thermostatted by the Nose-Hoover or Nose-Hoover chain method [P. Procacci et al., J. Chem. Phys. 125, 164101 (2006)]. That proof is essentially based on the fluctuation theorem by Evans and Searles [Adv. Phys. 51, 1529 (2002)] and on the equations of motion. Following an analogous approach, the authors derive here the Crooks equation in the context of molecular dynamics simulations of systems in the isothermal-isobaric (NPT) ensemble, whose dynamics is regulated by the Martyna-Tobias-Klein algorithm [J. Chem. Phys. 101, 4177 (1994)]. Their present derivation of the Crooks equation correlates to the demonstration of the Jarzynski identity for NPT systems recently proposed by Cuendet [J. Chem. Phys. 125, 144109 (2006)].     

Efficient use of nonequilibrium measurement to estimate free energy differences for molecular systems

A promising method for calculating free energy differences DeltaF is to generate nonequilibrium data via "fast-growth" simulations or by experiments--and then use Jarzynski's equality. However, a difficulty with using Jarzynski's equality is that DeltaF estimates converge very slowly and unreliably due to the nonlinear nature of the calculation--thus requiring large, costly data sets. The purpose of the work presented here is to determine the best estimate for DeltaF given a (finite) set of work values previously generated by simulation or experiment. Exploiting statistical properties of Jarzynski's equality, we present two fully automated analyses of nonequilibrium data from a toy model, and various simulated molecular systems. Both schemes remove at least several k(B)T of bias from DeltaF estimates, compared to direct application of Jarzynski's equality, for modest sized data sets (100 work values), in all tested systems. Results from one of the new methods suggest that good estimates of DeltaF can be obtained using 5-40-fold less data than was previously possible. Extending previous work, the new results exploit the systematic behavior of bias due to finite sample size. A key innovation is better use of the more statistically reliable information available from the raw data.     

Free energy calculation using molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. II. Thermodynamic integration along a spatial reaction coordinate

We propose the thermodynamic integration along a spatial reaction coordinate using the molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. This method provides a free energy calculation in solution along the reaction coordinate defined by the Cartesian coordinates of the solute atoms. The proposed method is based on the blue moon algorithm which can, in principle, handle any reaction coordinate as far as it is defined by the solute atom positions. In this article, we apply the present method to the complex formation process of the crown ether 18-Crown-6 (18C6) with the potassium ion in an aqueous solution. The separation between the geometric centers of these two molecules is taken to be the reaction coordinate for this system. The potential of mean force (PMF) becomes the maximum at the separation between the molecular centers being ～4 A?, which can be identified as the free energy barrier in the process of the molecular recognition. In a separation further than the free energy barrier, the PMF is slightly reduced to exhibit a plateau. In the region closer than the free energy barrier, approach of the potassium ion to the center of 18C6 also decreases the PMF. When the potassium ion is accommodated at the center of 18C6, the free energy is lower by -5.7 ± 0.7 kcal/mol than that at the above mentioned plateau or converged state. By comparing the results with those from the free energy calculation along the coupling parameters obtained in our previous paper [T. Miyata, Y. Ikuta, and F. Hirata, J. Chem. Phys. 133, 044114 (2010)], it is found that the effective interaction in water between 18C6 and the potassium ion vanishes beyond the molecular-center-separation of 10 A?. Furthermore, the conformation of 18C6 is found to be significantly changed depending upon the 18C6-K(+) distance. A proper conformational sampling and an accurate solvent treatment are crucial for realizing the accurate PMF, and we believe that the proposed method is useful to evaluate the PMF in a solution. A discussion upon the PMF in terms of the three-dimensional distribution function for the solvent is also presented.     

Absolute free energy calculations by thermodynamic integration in four spatial dimensions

An optimized technique for calculating the excess chemical potential of small molecules in dense liquids and the binding affinity of molecular ligands to biomolecules is reported. In this method, a molecular species is coupled to the system of interest via a nonphysical fourth spatial dimension w through which insertion or extraction can be carried out [R. Pomes, E. Eisenmesser, C. B. Post et al., J. Chem. Phys. 111, 3387 (1999)]. Molecular simulations are used to compute the potential of mean force (PMF) acting on the solute molecule in the fourth dimension. The excess chemical potential of that molecule is obtained as the difference in the PMF between fully coupled and fully decoupled systems. The simplicity, efficiency, and generality of the method are demonstrated for the calculation of the hydration free energies of water and methanol as well as sodium, cesium, and chloride ions. A significant advantage over other methods is that the 4D-PMF approach provides a single effective and general route for decoupling all nonbonded interactions (i.e., both Lennard-Jones and Coulombic) at once for both neutral and charged solutes. Direct calculation of the mean force from thermodynamic integration is shown to be more computationally efficient than calculating the PMF from umbrella sampling. Statistical error analysis suggests a simple strategy for optimizing sampling. The detailed analysis of systematic errors arising from the truncation of Coulombic interactions in a solvent droplet of finite size leads to straightforward corrections to ionic hydration free energies.     

The rectified second law of thermodynamics

Equilibrium thermodynamics is combined with Jarzynski's irreversible work theorem to quantify the excess entropy produced by irreversible processes. The resulting rectified form of the second law parallels the first law, in the sense that it facilitates the experimental measurement of excess entropy changes resulting from irreversible work and heat exchanges, just as the first law quantifies energy changes produced by either reversible or irreversible work and heat exchanges. The general form of the rectified second law is further applied to a broad class of quasi-static irreverisble (QSI) processes, for which all of the thermodynamic functions of both the system and surroundings remain continuously well-defined, thus facilitating excess entropy measurements by integrating exact differential functions along QSI paths. The results are illustrated by calculating the mechanical and thermal excess entropy produced by the irreversible unfolding of an RNA molecule.     

Activation energy for a model ferrous-ferric half reaction from transition path sampling

Activation parameters for the model oxidation half reaction of the classical aqueous ferrous ion are compared for different molecular simulation techniques. In particular, activation free energies are obtained from umbrella integration and Marcus theory based thermodynamic integration, which rely on the diabatic gap as the reaction coordinate. The latter method also assumes linear response, and both methods obtain the activation entropy and the activation energy from the temperature dependence of the activation free energy. In contrast, transition path sampling does not require knowledge of the reaction coordinate and directly yields the activation energy [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Benchmark activation energies from transition path sampling agree within statistical uncertainty with activation energies obtained from standard techniques requiring knowledge of the reaction coordinate. In addition, it is found that the activation energy for this model system is significantly smaller than the activation free energy for the Marcus model, approximately half the value, implying an equally large entropy contribution.     

Calculating potentials of mean force from steered molecular dynamics simulations

Steered molecular dynamics (SMD) permits efficient investigations of molecular processes by focusing on selected degrees of freedom. We explain how one can, in the framework of SMD, employ Jarzynski's equality (also known as the nonequilibrium work relation) to calculate potentials of mean force (PMF). We outline the theory that serves this purpose and connects nonequilibrium processes (such as SMD simulations) with equilibrium properties (such as the PMF). We review the derivation of Jarzynski's equality, generalize it to isobaric--isothermal processes, and discuss its implications in relation to the second law of thermodynamics and computer simulations. In the relevant regime of steering by means of stiff springs, we demonstrate that the work on the system is Gaussian-distributed regardless of the speed of the process simulated. In this case, the cumulant expansion of Jarzynski's equality can be safely terminated at second order. We illustrate the PMF calculation method for an exemplary simulation and demonstrate the Gaussian nature of the resulting work distribution.     

A water-swap reaction coordinate for the calculation of absolute protein-ligand binding free energies

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.     

Ordering of limits in the Jarzynski equality

We consider the sampling problems encountered in computing free-energy differences using Jarzynski's nonequilibrium work relation [Phys. Rev. Lett. 56, 2690 (1997)]. This relation expresses the free-energy change of a system, on which finite-time work is done, as an average over all possible trajectories of the system. This average can then be expressed as a cumulant expansion of the work. We study the scaling of these cumulants with an appropriately defined measure of phase-space accessibility epsilon and particle number N for two simple changes in state. We find that the cumulant expansion is slowly convergent for predominantly entropic processes, those where epsilon is considerably altered during the course of the process. An accurate determination of the free-energy change requires some knowledge of the behavior of the tails of the work distribution associated with the process. Jarzynski's irreversible work relation is only valid with the correct ordering of the infinite limits of N and epsilon, clarifying the regime of its applicability.     

Rovibrational spectroscopy calculations of neon dimer using a phase space truncated Weyl-Heisenberg wavelet basis

In a series of earlier articles [B. Poirier J. Theor. Comput. Chem. 2, 65 (2003); B. Poirier and A. Salam J. Chem. Phys. 121, 1690 (2004); B. Poirier and A. Salam J. Chem. Phys. 121, 1740 (2004)], a new method was introduced for performing exact quantum dynamics calculations in a manner that formally defeats exponential scaling with system dimensionality. The method combines an optimally localized, orthogonal Weyl-Heisenberg wavelet basis set with a simple phase space truncation scheme, and has already been applied to model systems up to 17 degrees of freedom (DOF's). In this paper, the approach is applied for the first time to a real molecular system (neon dimer), necessitating the development of an efficient numerical scheme for representing arbitrary potential energy functions in the wavelet representation. All bound rovibrational energy levels of neon dimer are computed, using both one DOF radial coordinate calculations and a three DOF Cartesian coordinate calculation. Even at such low dimensionalities, the approach is found to be competitive with another state-of-the-art method applied to the same system [J. Montgomery and B. Poirier J. Chem. Phys. 119, 6609 (2003)].     

Free energy landscapes of electron transfer system in dipolar environment below and above the rotational freezing temperature

Electron transfer reaction in a polar solvent is modeled by a solute dipole surrounded by dipolar molecules with simple rotational dynamics posted on the three-dimensional distorted lattice sites. The interaction energy between the solute and solvent dipoles as a reaction coordinate is adopted and free energy landscapes are calculated by generating all possible states for a 26 dipolar system and by employing Wang-Landau sampling algorithm for a 92 dipolar system. For temperatures higher than the energy scale of dipole-dipole interactions, the free energy landscapes for the small reaction coordinate region have quadratic shape as predicted by Marcus [Rev. Mod. Phys. 65, 599 (1993)] whereas for the large reaction coordinate region, the landscapes exhibit a nonquadratic shape. When the temperature drops, small notched structures appear on the free energy profiles because of the frustrated interactions among dipoles. The formation of notched structure is analyzed with statistical approach and it is shown that the amplitude of notched structure depend upon the segment size of the reaction coordinate and is characterized by the interaction energy among the dipoles. Using simulated free energy landscapes, the authors calculate the reaction rates as a function of the energy gap for various temperatures. At high temperature, the reactions rates follow a bell shaped (inverted parabolic) energy gap law in the small energy gap regions, while it becomes steeper than the parabolic shape in a large energy gap regions due to the nonquadratic shape of the free energy landscape. The peak position of parabola also changes as the function of temperature. At low temperature, the profile of the reaction rates is no longer smooth because of the many local minima of the free energy landscape.     

Adaptive steered molecular dynamics: validation of the selection criterion and benchmarking energetics in vacuum

The potential of mean force (PMF) for stretching decaalanine in vacuum was determined earlier by Park and Schulten [J. Chem. Phys. 120, 5946 (2004)] in a landmark article demonstrating the efficacy of combining steered molecular dynamics and Jarzynski's nonequilibrium relation. In this study, the recently developed adaptive steered molecular dynamics (ASMD) algorithm [G. Ozer, E. Valeev, S. Quirk, and R. Hernandez, J. Chem. Theory Comput. 6, 3026 (2010)] is used to reproduce the PMF of the unraveling of decaalanine in vacuum by averaging over fewer nonequilibrium trajectories. The efficiency and accuracy of the method are demonstrated through the agreement with the earlier work by Park and Schulten, a series of convergence checks compared to alternate SMD pulling strategies, and an analytical proof. The nonequilibrium trajectories obtained through ASMD have also been used to analyze the intrapeptide hydrogen bonds along the stretching coordinate. As the decaalanine helix is stretched, the initially stabilized i → i + 4 contacts (α-helix) is replaced by i → i + 3 contacts (3(10)-helix). No significant formation of i → i + 5 hydrogen bonds (π-helix) is observed.     

Potential of mean force calculations of ligand binding to ion channels from Jarzynski's equality and umbrella sampling

Potential of mean force (PMF) calculations provide a reliable method for determination of the absolute binding free energies for protein-ligand systems. The common method used for this purpose -- umbrella sampling with weighted histogram analysis -- is computationally very laborious, which limits its applications. Recently, a much simpler alternative for PMF calculations has become available, namely, using Jarzynski's equality in steered molecular dynamics simulations. So far, there have been a few comparisons of the two methods and mostly in simple systems that do not reflect the complexities of protein-ligand systems. Here, we use both methods to calculate the PMF for ion permeation and ligand binding to ion channels. Comparison of results indicate that Jarzynski's method suffers from relaxation problems in complex systems and would require much longer simulation times to yield reliable PMFs for protein-ligand systems.     

Springs and speeds in free energy reconstruction from irreversible single-molecule pulling experiments

The nonequilibrium work relation allows for the calculation of equilibrium free energy differences between states based on the exponential average of accumulated work from irreversible transitions. Here, we compare two distinct approaches of calculating free energy surfaces from unidirectional single-molecule pulling experiments: the stiff spring approximation and the Hummer-Szabo method. First, we perform steered molecular dynamics simulations to mechanically stretch the model peptide deca-alanine using harmonic potentials with different spring stiffnesses and at various constant pulling velocities. We then calculate free energy surfaces based on the two methods and their variants, including the first and second cumulant expansion of the exponentially weighted work and the Gaussian position approximation for the delta function in Hummer and Szabo's expression. We find that with large harmonic force constants, the second cumulant expansion performs well in conjunction with either the stiff spring approximation or the Hummer-Szabo method. When interpreting dynamic force spectroscopy (pullings at different speeds), the second cumulant expansion of the stiff spring approximation performs the best when pulling velocities are similar, but variants of the Hummer-Szabo perform the best when they are spread over a large spectrum. While these conclusion are not definitive for all systems, the insights should prove useful for scientists interpreting nonequilibrium pulling experiments.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode

The thermochemical constants for the oxidation of tyrosine and tryptophan through proton coupled electron transfer in aqueous solution have been computed applying a recently developed density functional theory (DFT) based molecular dynamics method for reversible elimination of protons and electrons. This method enables us to estimate the solvation free energy of a proton (H(+)) in a periodic model system from the free energy for the deprotonation of an aqueous hydronium ion (H(3)O(+)). Using the computed solvation free energy of H(+) as reference, the deprotonation and oxidation free energies of an aqueous species can be converted to pK(a) and normal hydrogen electrode (NHE) potentials. This conversion requires certain thermochemical corrections which were first presented in a similar study of the oxidation of hydrobenzoquinone [J. Cheng, M. Sulpizi, and M. Sprik, J. Chem. Phys. 131, 154504 (2009)]. Taking a different view of the thermodynamic status of the hydronium ion, these thermochemical corrections are revised in the present work. The key difference with the previous scheme is that the hydronium is now treated as an intermediate in the transfer of the proton from solution to the gas-phase. The accuracy of the method is assessed by a detailed comparison of the computed pK(a), NHE potentials and dehydrogenation free energies to experiment. As a further application of the technique, we have analyzed the role of the solvent in the oxidation of tyrosine by the tryptophan radical. The free energy change computed for this hydrogen atom transfer reaction is very similar to the gas-phase value, in agreement with experiment. The molecular dynamics results however, show that the minimal solvent effect on the reaction free energy is accompanied by a significant reorganization of the solvent.     

Potentials of mean force for acetylcholine unbinding from the alpha7 nicotinic acetylcholine receptor ligand-binding domain

The nicotinic acetylcholine receptor is a prototype ligand-gated ion channel that mediates signal transduction in the neuromuscular junctions and other cholinergic synapses. The molecular basis for the energetics of ligand binding and unbinding is critical to our understanding of the pharmacology of this class of receptors. Here, we used steered molecular dynamics to investigate the unbinding of acetylcholine from the ligand-binding domain of human alpha7 nicotinic acetylcholine receptor along four different predetermined pathways. Pulling forces were found to correlate well with interactions between acetylcholine and residues in the binding site during the unbinding process. From multiple trajectories along these unbinding pathways, we calculated the potentials of mean force for acetylcholine unbinding. Four available methods based on Jarzynski's equality were used and compared for their efficiencies. The most probable pathway was identified to be along a direction approximately parallel to the membrane. The derived binding energy for acetylcholine was in good agreement with that derived from the experimental binding constant for acetylcholine binding protein, but significantly higher than that for the complete human alpha7 nicotinic acetylcholine receptor. In addition, it is likely that several intermediate states exist along the unbinding pathways.