A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions

The recently introduced hills method (Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562) is a powerful tool to compute the multidimensional free energy surface of intrinsically concerted reactions. We have extended this method by focusing our attention on localizing the lowest free energy path that connects the stable reactant and product states. This path represents the most probable reaction mechanism, similar to the zero temperature intrinsic reaction coordinate, but also includes finite temperature effects. The transformation of the multidimensional problem to a one-dimensional reaction coordinate allows for accurate convergence of the free energy profile along the lowest free energy path using standard free energy methods. Here we apply the hills method, our lowest free energy path search algorithm, and umbrella sampling to the prototype S(N)2 reaction. The hills method replaces the in many cases difficult problem of finding a good reaction coordinate with choosing relatively simple collective variables, such as the bond lengths of the broken and formed chemical bonds. The second part of the paper presents a guide to using the hills method, in which we test and fine-tune the method for optimal accuracy and efficiency using the umbrella sampling results as a reference.     

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.     

The thermal stability of lattice-energy minima of 5-fluorouracil: metadynamics as an aid to polymorph prediction

This paper reports a novel methodology for the free-energy minimization of crystal structures exhibiting strong, anisotropic interactions due to hydrogen bonding. The geometry of the thermally expanded cell was calculated by exploiting the dependence of the free-energy derivatives with respect to cell lengths and angles on the average pressure tensor computed in short molecular dynamics simulations. All dynamic simulations were performed with an elaborate anisotropic potential based on a distributed multipole analysis of the isolated molecule charge density. Changes in structure were monitored via simulated X-ray diffraction patterns. The methodology was used to minimize the free energy at ambient conditions of a set of experimental and hypothetical 5-fluorouracil crystal structures, generated in a search for lattice-energy minima with the same model potential. Our results demonstrate that the majority ( approximately 75%) of lattice-energy minima are thermally stable at ambient conditions, and hence, the free-energy (like the lattice-energy) surface is complex and highly undulating. Metadynamics trajectories (Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562) started from the free-energy minima only produced transitions that preserved the hydrogen-bonding motif, and thus, further developments are needed for this method to efficiently explore such free-energy surfaces. The existence of so many free-energy minima, with large barriers for the alteration of the hydrogen-bonding motif, is consistent with the range of motifs observed in crystal structures of 5-fluorouracil and other 5-substituted uracils.     

Flexible docking in solution using metadynamics

We apply our recently developed metadynamics method to the docking of ligands on flexible receptors in water solution. This method mimics the real dynamics of a ligand exiting or entering an enzyme and in so doing reconstructs the free energy surface. We apply it to four docking cases: beta-trypsin/benzamidine, beta-trypsin/chlorobenzamidine, immunoglobulin McPC-603/phosphocholine, and cyclin-dependent kinase 2/staurosporine. In every case studied, the method is able to predict the docked geometry and the free energy of docking. Its added value with respect to many other available methods is that it reconstructs the complete free energy surface, including all the relevant minima and the barriers between them.     

Computing accurate potentials of mean force in electrolyte solutions with the generalized gradient-augmented harmonic Fourier beads method

We establish the accuracy of the novel generalized gradient-augmented harmonic Fourier beads (ggaHFB) method in computing free-energy profiles or potentials of mean force (PMFs) through comparison with two independent conventional techniques. In particular, we employ umbrella sampling with one dimensional weighted histogram analysis method (WHAM) and free molecular dynamics simulation of radial distribution functions to compute the PMF for the Na(+)-Cl(-) ion-pair separation to 16 A in 1.0M NaCl solution in water. The corresponding ggaHFB free-energy profile in six dimensional Cartesian space is in excellent agreement with the conventional benchmarks. We then explore changes in the PMF in response to lowering the NaCl concentration to physiological 0.3 and 0.1M, and dilute 0.0M concentrations. Finally, to expand the scope of the ggaHFB method, we formally develop the free-energy gradient approximation in arbitrary nonlinear coordinates. This formal development underscores the importance of the logarithmic Jacobian correction to reconstruct true PMFs from umbrella sampling simulations with either WHAM or ggaHFB techniques when nonlinear coordinate restraints are used with Cartesian propagators. The ability to employ nonlinear coordinates and high accuracy of the computed free-energy profiles further advocate the use of the ggaHFB method in studies of rare events in complex systems.     

Metadynamics in essential coordinates: free energy simulation of conformational changes

We propose an approach that combines an extraction of collective motions of a molecular system with a sampling of its free energy surface. A recently introduced method of metadynamics allows exploration of the free energy surface of a molecular system by means of coarse-grained dynamics with flooding of free energy minima. This free energy surface is defined as a function of a set of collective variables (e.g., interatomic distances, angles, torsions, and others). In this study, essential coordinates determined by essential dynamics (principle component analysis) were used as collective variables in metadynamics. First, dynamics of the model system (explicitly solvated alanine dipeptide, Ace-Ala-Nme) was simulated by a classical molecular dynamics simulation. The trajectory (1 ns) was then analyzed by essential dynamics to obtain essential coordinates. The free energy surface as a function of the first and second essential coordinates was then explored by metadynamics. The resulting free energy surface is in agreement with other studies of this system. We propose that a combination of these two methods (metadynamics and essential dynamics) has great potential in studies of conformational changes in peptides and proteins.     

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.     

Harmonic Fourier beads method for studying rare events on rugged energy surfaces

We present a robust, distributable method for computing minimum free energy paths of large molecular systems with rugged energy landscapes. The method, which we call harmonic Fourier beads (HFB), exploits the Fourier representation of a path in an appropriate coordinate space and proceeds iteratively by evolving a discrete set of harmonically restrained path points-beads-to generate positions for the next path. The HFB method does not require explicit knowledge of the free energy to locate the path. To compute the free energy profile along the final path we employ an umbrella sampling method in two generalized dimensions. The proposed HFB method is anticipated to aid the study of rare events in biomolecular systems. Its utility is demonstrated with an application to conformational isomerization of the alanine dipeptide in gas phase.     

Tracing the minimum-energy path on the free-energy surface

The free-energy profile of a reaction can be estimated in a molecular-dynamics approach by imposing a mechanical constraint along a reaction coordinate (RC). Many recent studies have shown that the temperature can greatly influence the path followed by the reactants. Here, we propose a practical way to construct the minimum-energy path directly on the free-energy surface at a given temperature. First, we follow the blue-moon ensemble method to derive the expression of the free-energy gradient for a given RC. These derivatives are then used to find the actual minimum-energy reaction path at finite temperature, in a way similar to the intrinsic reaction path of Fukui on the potential-energy surface. [K. Fukui, J. Phys. Chem. 74, 4161 (1970)]. Once the path is known, one can calculate the free-energy profile using thermodynamic integration. We also show that the mass-metric correction cancels for many types of constraints, making the procedure easy to use. Finally, the minimum-free-energy path at 300 K for the addition of CCl2 to ethylene is compared with a path based on a simple one-dimensional reaction coordinate. A comparison is also given with the reaction path at 0 K.     

Umbrella Sampling Simulations of Carbon Nanoparticles Crossing Immiscible Solvents

We use molecular dynamics to compute the free energy of carbon nanoparticles crossing a hydrophobic–hydrophilic interface. The simulations are performed on a biphasic system consisting of immiscible solvents (i.e., cyclohexane and water). We solvate a carbon nanoparticle into the cyclohexane layer and use a pull force to drive the nanoparticle into water, passing over the interface. Next, we accumulate a series of umbrella sampling simulations along the path of the nanoparticle and compute the solvation free energy with respect to the two solvents. We apply the method on three carbon nanoparticles (i.e., a carbon nanocone, a nanotube, and a graphene nanosheet). In addition, we record the water-accessible surface area of the nanoparticles during the umbrella simulations. Although we detect complete wetting of the external surface of the nanoparticles, the internal surface of the nanotube becomes partially wet, whereas that of the nanocone remains dry. This is due to the nanoconfinement of the particular nanoparticles, which shields the hydrophobic interactions encountered inside the pores. We show that cyclohexane molecules remain attached on the concave surface of the nanotube or the nanocone without being disturbed by the water molecules entering the cavity.     

From A to B in free energy space

The authors present a new method for searching low free energy paths in complex molecular systems at finite temperature. They introduce two variables that are able to describe the position of a point in configurational space relative to a preassigned path. With the help of these two variables the authors combine features of approaches such as metadynamics or umbrella sampling with those of path based methods. This allows global searches in the space of paths to be performed and a new variational principle for the determination of low free energy paths to be established. Contrary to metadynamics or umbrella sampling the path can be described by an arbitrary large number of variables, still the energy profile along the path can be calculated. The authors exemplify the method numerically by studying the conformational changes of alanine dipeptide.     

Automated Search of Minimum Free‐Energy Path by Umbrella Integration

The free‐energy landscape is an important factor for understanding the conformational equilibria of chemical reactions, and many techniques have been developed to calculate the potential of the mean force. Unfortunately, these methods require a previous knowledge of the system for calculations because the results depend on the reaction coordinates. In this study, we combine the scaled hypersphere search method with the umbrella integration method to obtain the transition states on free‐energy landscapes and minimum‐free‐energy paths (MFEPs). With this approach, the MFEP connections between known and unknown equilibrium points are constructed without the prior knowledge of the free‐energy landscape. The problem of reaction coordinates can be solved by using a multidimensional, fully automated interrogation of MFEPs for acquiring the potential of mean force. The efficiency of the proposed method is demonstrated by applying it to alanine dipeptide and alanine tripeptide. © 2018 Wiley Periodicals, Inc.     

Structural and pathway complexity of beta-strand reorganization within aggregates of human transthyretin(105-115) peptide

Interstrand conformational rearrangements of human transthyretin peptide (TTR(105-115)) within dimeric aggregates were simulated by means of molecular dynamics (MD) with implicit solvation model for a total length of 48 micros. The conformations sampled in the MD simulations were clustered to identify free energy minima without any projections of free energy surface. A connected graph was constructed with nodes (=clusters) and edges corresponding to free energy minima and transitions between nodes, respectively. This connected graph which reflects the complexity of the free energy surface was used to extract the transition disconnectivity graph, which reflects the overall free energy barriers between pairs of free energy minima but does not contain information on transition paths. The routes of transitions between important free energy minima were obtained by further processing the original graph and the MD data. We have found that both parallel and antiparallel aggregates are populated. The parallel aggregates with different alignment patterns are separated by nonnegligible free energy barriers. Multiroutes exist in the interstrand conformational reorganization. Most visited routes do not dominant the kinetics, while less visited routes contribute a little each but they are numerous and their total contributions are actually dominant. There are various kinds of reptation motions, including those through a beta-bulge, side-chain aided reptation, and flipping or rotation of a hairpin formed by one strand.     

Protein structure prediction using basin-hopping

Associative memory Hamiltonian structure prediction potentials are not overly rugged, thereby suggesting their landscapes are like those of actual proteins. In the present contribution we show how basin-hopping global optimization can identify low-lying minima for the corresponding mildly frustrated energy landscapes. For small systems the basin-hopping algorithm succeeds in locating both lower minima and conformations closer to the experimental structure than does molecular dynamics with simulated annealing. For large systems the efficiency of basin-hopping decreases for our initial implementation, where the steps consist of random perturbations to the Cartesian coordinates. We implemented umbrella sampling using basin-hopping to further confirm when the global minima are reached. We have also improved the energy surface by employing bioinformatic techniques for reducing the roughness or variance of the energy surface. Finally, the basin-hopping calculations have guided improvements in the excluded volume of the Hamiltonian, producing better structures. These results suggest a novel and transferable optimization scheme for future energy function development.     

On the efficient path integral evaluation of thermal rate constants within the quantum instanton approximation

We present an efficient path integral approach for evaluating thermal rate constants within the quantum instanton (QI) approximation that was recently introduced to overcome the quantitative deficiencies of the earlier semiclassical instanton approach [Miller, Zhao, Ceotto, and Yang, J. Chem. Phys. 119, 1329 (2003)]. Since the QI rate constant is determined solely by properties of the (quantum) Boltzmann operator (specifically, by the zero time properties of the flux-flux and delta-delta correlation functions), it can be evaluated by well-established techniques of imaginary time path integrals even for quite complex chemical reactions. Here we present a series of statistical estimators for relevant quantities which can be evaluated straightforwardly with any nonlinear reaction coordinates and general Hamiltonians in Cartesian space. To facilitate the search for the optimal dividing surfaces required by the QI approximation, we introduce a two-dimensional quantum free energy surface associated with the delta-delta correlation function and describe how an adaptive umbrella sampling can be used effectively to construct such a free energy surface. The overall computational procedure is illustrated by the application to a hydrogen exchange reaction in gas phase, which shows excellent agreement of the QI rates with those obtained from quantum scattering calculations.     

Convergence of Free Energy Profile of Coumarin in Lipid Bilayer

Atomistic molecular dynamics (MD) simulations of druglike molecules embedded in lipid bilayers are of considerable interest as models for drug penetration and positioning in biological membranes. Here we analyze partitioning of coumarin in dioleoylphosphatidylcholine (DOPC) bilayer, based on both multiple, unbiased 3 μs MD simulations (total length) and free energy profiles along the bilayer normal calculated by biased MD simulations (～7 μs in total). The convergences in time of free energy profiles calculated by both umbrella sampling and z-constraint techniques are thoroughly analyzed. Two sets of starting structures are also considered, one from unbiased MD simulation and the other from "pulling" coumarin along the bilayer normal. The structures obtained by pulling simulation contain water defects on the lipid bilayer surface, while those acquired from unbiased simulation have no membrane defects. The free energy profiles converge more rapidly when starting frames from unbiased simulations are used. In addition, z-constraint simulation leads to more rapid convergence than umbrella sampling, due to quicker relaxation of membrane defects. Furthermore, we show that the choice of RESP, PRODRG, or Mulliken charges considerably affects the resulting free energy profile of our model drug along the bilayer normal. We recommend using z-constraint biased MD simulations based on starting geometries acquired from unbiased MD simulations for efficient calculation of convergent free energy profiles of druglike molecules along bilayer normals. The calculation of free energy profile should start with an unbiased simulation, though the polar molecules might need a slow pulling afterward. Results obtained with the recommended simulation protocol agree well with available experimental data for two coumarin derivatives.     

Multicoordinate driven method for approximating enzymatic reaction paths: automatic definition of the reaction coordinate using a subset of chemical coordinates

We present a generalization of the reaction coordinate driven method to find reaction paths and transition states for complicated chemical processes, especially enzymatic reactions. The method is based on the definition of a subset of chemical coordinates; it is simple, robust, and suitable to calculate one or more alternative pathways, intermediate minima, and transition-state geometries. Though the results are approximate and the computational cost is relatively high, the method works for large systems, where others often fail. It also works when a certain reaction path competes with others having a lower energy barrier. Accordingly, the procedure is appropriate to test hypothetical reaction mechanisms for complicated systems and provides good initial guesses for more accurate methods. We present tests on a number of simple reactions and on several complicated chemical transformations and compare the results with those obtained by other methods. Calculation of the reaction path for the enzymatic hydrolysis of the substrate by dUTPase for an active-site model with 85 atoms, including several loosely bound water molecules, indicates that the method is feasible for the study of enzyme mechanisms.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

Free energy landscape of A-DNA to B-DNA conversion in aqueous solution

The interconversion between the well-characterized A- and B-forms of DNA is a structural transition for which the intermediate states and the free energy difference between the two endpoints are not known precisely. In the present study, the difference between the Root Mean Square Distance (RMSD) from canonical A-form and B-form DNA is used as an order parameter to characterize this free energy difference using umbrella sampling molecular dynamics (MD) simulations with explicit solvent. The constraint imposed along this order parameter allows relatively unrestricted evolution of the intermediate structures away from both canonical A- and B-forms. The free energy difference between the A- and B-forms for the hexamer DNA sequence CTCGAG in aqueous solution is conservatively estimated to be at least 2.8 kcal/mol. A continuum of intermediate structures with no well-defined local minima links the two forms. The absence of any major barriers in the free energy surface is consistent with spontaneous conversion of the A-form DNA to B-form DNA in unconstrained simulations. The extensive sampling in the MD simulations (>0.1 mus) also allowed quantitative energetic characterization of local backbone conformational variables such as sugar pseudorotation angles and BI/BII state equilibria and their dependence on base identity. The absolute minimum in the calculated free energy profile corresponds closely to the crystal structure of the hexamer sequence, indicating that the present method has the potential to identify the most stable state for an arbitrary DNA sequence in water.