Calculation of free energy landscapes: A histogram reweighted metadynamics approach

We present an efficient method for the calculation of free energy landscapes. Our approach involves a history‐dependent bias potential, which is evaluated on a grid. The corresponding free energy landscape is constructed via a histogram reweighting procedure a posteriori. Because of the presence of the bias potential, it can be also used to accelerate rare events. In addition, the calculated free energy landscape is not restricted to the actual choice of collective variables and can in principle be extended to auxiliary variables of interest without further numerical effort. The applicability is shown for several examples. We present numerical results for the alanine dipeptide and the Met‐Enkephalin in explicit solution to illustrate our approach. Furthermore, we derive an empirical formula that allows the prediction of the computational cost for the ordinary metadynamics variant in comparison with our approach, which is validated by a dimensionless representation. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Finite-Temperature Dimer Method for Finding Saddle Points on Free Energy Surfaces

The dimer method and its variants have been shown to be efficient in finding saddle points on potential surfaces. In the dimer method, the most unstable direction is approximately obtained by minimizing the total potential energy of the dimer. Then, the force in this direction is reversed to move the dimer toward saddle points. When the finite-temperature effect is important for a high-dimensional system, one usually needs to describe the dynamics in a low-dimensional space of reaction coordinates. In this case, transition states are collected as saddle points on the free energy surface. The traditional dimer method cannot be directly employed to find saddle points on a free energy surface since the surface is not known a priori. Here, we develop a finite-temperature dimer method for searching saddle points on the free energy surface. In this method, a constrained rotation dynamics of the dimer system is used to sample dimer directions and an efficient average method is used to obtain a good approximation of the most unstable direction. This approximated direction is then used in reversing the force component and evolving the dimer toward saddle points. Our numerical results suggest that the new method is efficient in finding saddle points on free energy surfaces. © 2019 Wiley Periodicals, Inc.     

Exploring reaction pathways with transition path and umbrella sampling: application to methyl maltoside

The transition path sampling (TPS) method is a powerful approach to study chemical reactions or transitional properties on complex potential energy landscapes. One of the main advantages of the method over potential of mean force methods is that reaction rates can be directly accessed without knowledge of the exact reaction coordinate. We have investigated the complementary nature of these two differing approaches, comparing transition path sampling with the weighted histogram analysis method to study a conformational change in a small model system. In this case study, the transition paths for a transition between two rotational conformers of a model disaccharide molecule, methyl beta-D-maltoside, were compared with a free energy surface constrained by the two commonly used glycosidic (phi,psi) torsional angles. The TPS method revealed a reaction channel that was not apparent from the potential of mean force method, and the suitability of phi and psi as reaction coordinates to describe the isomerization in vacuo was confirmed by examination of the transition path ensemble. Using both transition state theory and transition path sampling methods, the transition rate was estimated. We have estimated a characteristic time between transitions of approximately 160 ns for this rare isomerization event between the two conformations of the carbohydrate. We conclude that transition path sampling can extract subtle information about the dynamics not apparent from the potential of mean force method. However, in calculating the reaction rate, the transition path sampling method required 27.5 times the computational effort than was needed by the potential of mean force method.     

Implementation of an adaptive umbrella sampling method for the calculation of multidimensional potential of mean force of chemical reactions in solution

We describe the implementation of an adaptive umbrella sampling method, making use of the weighted histogram analysis method, for computing multidimensional potential of mean force for chemical reaction in solution. The approach is illustrated by investigating the effect of aqueous solution on the free energy surface for the proton transfer reaction of [H(3)N-H-NH(3)](+) using a combined quantum mechanical and molecular mechanical AM1/TIP3P potential.     

Efficient calculation of relative binding free energies by umbrella sampling perturbation

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide‐protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1–2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation. © 2014 Wiley Periodicals, Inc.     

Exploring chemical reactivity of complex systems with path‐based coordinates: Role of the distance metric

Path‐based reaction coordinates constitute a valuable tool for free‐energy calculations in complex processes. When a reference path is defined by means of collective variables, a nonconstant distance metric that incorporates the nonorthonormality of these variables should be taken into account. In this work, we show that, accounting for the correct metric tensor, these kind of variables can provide iso‐hypersurfaces that coincide with the iso‐committor surfaces and that activation free energies equal the value that would be obtained if the committor function itself were used as reaction coordinate. The advantages of the incorporation of the variable metric tensor are illustrated with the analysis of the enzymatic reaction catalyzed by isochorismate‐pyruvate lyase. Hybrid QM/MM techniques are used to obtain the free energy profile and to analyze reactive trajectories initiated at the transition state. For this example, the committor histogram is peaked at 0.5 only when a variable metric tensor is incorporated in the definition of the path‐based coordinate. © 2014 Wiley Periodicals, Inc.     

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.     

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‐fly reconstruction of free‐energy profiles using logarithmic mean‐force dynamics

Mean‐force dynamics (MFD), which is a fictitious dynamics for a set of collective variables on a potential of mean‐force, is a powerful algorithm to efficiently explore free‐energy landscapes. Recently, we have introduced logarithmic MFD (LogMFD) (Morishita et al., Phys. Rev. E 2012, 85, 066702) which overcomes difficulties encounterd in free‐energy calculations using standard approaches such as thermodynamic integration. Here, we present a guide to implementing LogMFD calculations paying attention to the practical issues in choosing the parameters in LogMFD. A primary focus is given to the effect of the parameters on the accuracy of the reconstructed free‐energy profiles. A recipe for reducing the errors due to energy dissipation is presented. We also demonstrate that multidimensional free‐energy landscapes can be reconstructed on‐the‐fly using LogMFD, which cannot be accomplished using any other free‐energy calculation techniques. © 2013 Wiley Periodicals, Inc.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 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.     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

Tinker‐OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs

The capabilities of the polarizable force fields for alchemical free energy calculations have been limited by the high computational cost and complexity of the underlying potential energy functions. In this work, we present a GPU‐based general alchemical free energy simulation platform for polarizable potential AMOEBA. Tinker‐OpenMM, the OpenMM implementation of the AMOEBA simulation engine has been modified to enable both absolute and relative alchemical simulations on GPUs, which leads to a ∼200‐fold improvement in simulation speed over a single CPU core. We show that free energy values calculated using this platform agree with the results of Tinker simulations for the hydration of organic compounds and binding of host–guest systems within the statistical errors. In addition to absolute binding, we designed a relative alchemical approach for computing relative binding affinities of ligands to the same host, where a special path was applied to avoid numerical instability due to polarization between the different ligands that bind to the same site. This scheme is general and does not require ligands to have similar scaffolds. We show that relative hydration and binding free energy calculated using this approach match those computed from the absolute free energy approach. © 2017 Wiley Periodicals, Inc.     

Single-sweep methods for free energy calculations

A simple, efficient, and accurate method is proposed to map multidimensional free energy landscapes. The method combines the temperature-accelerated molecular dynamics (TAMD) proposed in [L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006)] with a variational reconstruction method using radial-basis functions for the representation of the free energy. TAMD is used to rapidly sweep through the important regions of the free energy landscape and to compute the gradient of the free energy locally at points in these regions. The variational method is then used to reconstruct the free energy globally from the mean force at these points. The algorithmic aspects of the single-sweep method are explained in detail, and the method is tested on simple examples and used to compute the free energy of the solvated alanine dipeptide in two and four dihedral angles.     

Determination of conformational free energies of peptides by multidimensional adaptive umbrella sampling

We improve the multidimensional adaptive umbrella sampling method for the computation of conformational free energies of biomolecules. The conformational transition between the alpha-helical and beta-hairpin conformational states of an alanine decapeptide is used as an example. Convergence properties of the weighted-histogram-analysis-based adaptive umbrella sampling can be improved by using multiple replicas in each adaptive iteration and by using adaptive updating of the bounds of the umbrella potential. Using positional root-mean-square deviations from structures of the alpha-helical and beta-hairpin reference states as reaction coordinates, we obtained well-converged free energy surfaces of both the in-vacuum and in-solution decapeptide systems. From the free energy surfaces well-converged relative free energies between the two conformational states can be derived. Advantages and disadvantages of different methods for obtaining conformational free energies as well as implications of our results in studying conformational transitions of proteins and in improving force field are discussed.     

Virtual‐system‐coupled adaptive umbrella sampling to compute free‐energy landscape for flexible molecular docking

A novel enhanced conformational sampling method, virtual‐system‐coupled adaptive umbrella sampling (V‐AUS), was proposed to compute 300‐K free‐energy landscape for flexible molecular docking, where a virtual degrees of freedom was introduced to control the sampling. This degree of freedom interacts with the biomolecular system. V‐AUS was applied to complex formation of two disordered amyloid‐β (Aβ_30–35) peptides in a periodic box filled by an explicit solvent. An interpeptide distance was defined as the reaction coordinate, along which sampling was enhanced. A uniform conformational distribution was obtained covering a wide interpeptide distance ranging from the bound to unbound states. The 300‐K free‐energy landscape was characterized by thermodynamically stable basins of antiparallel and parallel β‐sheet complexes and some other complex forms. Helices were frequently observed, when the two peptides contacted loosely or fluctuated freely without interpeptide contacts. We observed that V‐AUS converged to uniform distribution more effectively than conventional AUS sampling did. © 2015 Wiley Periodicals, Inc.     

Producing DFT/MM enzyme reaction trajectories from SCC‐DFTB/MM driving forces to probe the underlying electronics of a glycosyltransferase reaction

The SCC‐DFTB/MIO/CHARMM free energy surface for a glycosyltransferase, TcTS, is benchmarked against a DFT/MM reaction trajectory using the same CHARMM MM force field ported to the NWChem package. The popular B3LYP functional, against which the MIO parameter set was parameterized is used to optimize TS structures and run DFT reaction dynamics. A novel approach was used to generate reaction forces from a SCC‐DFTB/MIO/CHARMM reaction surface to drive B3LYP/6‐31G/MM and B3LYP/6‐31G(d)/MM reaction trajectories. Although TS structures compare favorably, differences stemming primarily from a minimal basis set approximation prevented a successful 6‐31G(d) FEARCF reaction dynamics trajectory. None the less, the dynamic evolution of the B3LYP/6‐31G/MM‐computed electron density provided an opportunity to perform NBO analysis along the reaction trajectory. Here, we illustrate that a successful ab initio reaction trajectory is computationally accessible when the underlying potential energy function of the semi‐empirical method used to produce driving forces is sufficiently close to the ab initio potential. © 2017 Wiley Periodicals, Inc.     

Topography of the free-energy landscape probed via mechanical unfolding of proteins

Single-molecule experiments in which proteins are unfolded by applying mechanical stretching forces generally force unfolding to proceed along a reaction coordinate that is different from that in chemical or thermal denaturation. Here we simulate the mechanical unfolding and refolding of a minimalist off-lattice model of the protein ubiquitin to explore in detail the slice of the multidimensional free-energy landscape that is accessible via mechanical pulling experiments. We find that while the free-energy profile along typical "chemical" reaction coordinates may exhibit two minima, corresponding to the native and denatured states, the free energy G(z) is typically a monotonic function of the mechanical coordinate z equal to the protein extension. Application of a stretching force along z tilts the free-energy landscape resulting in a bistable (or multistable) free energy G(z)-fz probed in mechanical unfolding experiments. We construct a two-dimensional free-energy surface as a function of both chemical and mechanical reaction coordinates and examine the coupling between the two. We further study the refolding trajectories after the protein has been prestretched by a large force, as well as the mechanical unfolding trajectories in the presence of a large stretching force. We demonstrate that the stretching forces required to destabilize the native state thermodynamically are larger than those expected on the basis of previous experimental estimates of G(z). This finding is consistent with the recent experimental studies, indicating that proteins may refold even in the presence of a substantial stretching force. Finally, we show that for certain temperatures the free energy of a polyprotein chain consisting of multiple domains is a linear function of the chain extension. We propose that the recently observed "slow phase" in the refolding of proteins under mechanical tension may be viewed as downhill diffusion in such a linear potential.     

Absolute free energies of biomolecules from unperturbed ensembles

Computing the absolute free energy of a macromolecule's structural state, F, is a challenging problem of high relevance. This study presents a method that computes F using only information from an unperturbed simulation of the macromolecule in the relevant conformational state, ensemble, and environment. Absolute free energies produced by this method, dubbed V aluation of L ocal C onfiguration I ntegral with D ynamics (VALOCIDY), enable comparison of alternative states. For example, comparing explicitly solvated and vaporous states of amino acid side‐chain analogs produces solvation free energies in good agreement with experiments. Also, comparisons between alternative conformational states of model heptapeptides (including the unfolded state) produce free energy differences in agreement with data from μs molecular‐dynamics simulations and experimental propensities. The potential of using VALOCIDY in computational protein design is explored via a small design problem of stabilizing a β‐turn structure. When VALOCIDY‐based estimation of folding free energy is used as the design metric, the resulting sequence folds into the desired structure within the atomistic force field used in design. The VALOCIDY‐based approach also recognizes the distinct status of the native sequence regardless of minor details of the starting template structure, in stark contrast with a traditional fixed‐backbone approach. © 2013 Wiley Periodicals, Inc.