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.     

Adaptive‐numerical‐bias metadynamics

A metadynamics scheme is presented in which the free energy surface is filled with progressively adding adaptive biasing potentials, obtained from the accumulated probability distribution of the collective variables. Instead of adding Gaussians with assigned height and width in conventional metadynamics method, here we add a more realistic adaptive biasing potential to the Hamiltonian of the system. The shape of the adaptive biasing potential is adjusted on the fly by sampling over the visited states. As the top of the barrier is approached, the biasing potentials become wider. This decreases the problem of trapping the system in the niches, introduced by the addition of Gaussians of fixed height in metadynamics. Our results for the free energy profiles of three test systems show that this method is more accurate and converges more quickly than the conventional metadynamics, and is quite comparable (in accuracy and convergence rate) with the well‐tempered metadynamics method. © 2017 Wiley Periodicals, Inc.     

Sampling free energy surfaces as slices by combining umbrella sampling and metadynamics

Metadynamics (MTD) is a very powerful technique to sample high‐dimensional free energy landscapes, and due to its self‐guiding property, the method has been successful in studying complex reactions and conformational changes. MTD sampling is based on filling the free energy basins by biasing potentials and thus for cases with flat, broad, and unbound free energy wells, the computational time to sample them becomes very large. To alleviate this problem, we combine the standard Umbrella Sampling (US) technique with MTD to sample orthogonal collective variables (CVs) in a simultaneous way. Within this scheme, we construct the equilibrium distribution of CVs from biased distributions obtained from independent MTD simulations with umbrella potentials. Reweighting is carried out by a procedure that combines US reweighting and Tiwary–Parrinello MTD reweighting within the Weighted Histogram Analysis Method (WHAM). The approach is ideal for a controlled sampling of a CV in a MTD simulation, making it computationally efficient in sampling flat, broad, and unbound free energy surfaces. This technique also allows for a distributed sampling of a high‐dimensional free energy surface, further increasing the computational efficiency in sampling. We demonstrate the application of this technique in sampling high‐dimensional surface for various chemical reactions using ab initio and QM/MM hybrid molecular dynamics simulations. Further, to carry out MTD bias reweighting for computing forward reaction barriers in ab initio or QM/MM simulations, we propose a computationally affordable approach that does not require recrossing trajectories. © 2016 Wiley Periodicals, Inc.     

BAR-based optimum adaptive steered MD for configurational sampling

The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.     

Using the local elevation method to construct optimized umbrella sampling potentials: Calculation of the relative free energies and interconversion barriers of glucopyranose ring conformers in water

A method is proposed to combine the local elevation (LE) conformational searching and the umbrella sampling (US) conformational sampling approaches into a single local elevation umbrella sampling (LEUS) scheme for (explicit‐solvent) molecular dynamics (MD) simulations. In this approach, an initial (relatively short) LE build‐up (searching) phase is used to construct an optimized biasing potential within a subspace of conformationally relevant degrees of freedom, that is then used in a (comparatively longer) US sampling phase. This scheme dramatically enhances (in comparison with plain MD) the sampling power of MD simulations, taking advantage of the fact that the preoptimized biasing potential represents a reasonable approximation to the negative of the free energy surface in the considered conformational subspace. The method is applied to the calculation of the relative free energies of β‐D ‐glucopyranose ring conformers in water (within the GROMOS 45A4 force field). Different schemes to assign sampled conformational regions to distinct states are also compared. This approach, which bears some analogies with adaptive umbrella sampling and metadynamics (but within a very distinct implementation), is shown to be: (i) efficient (nearly all the computational effort is invested in the actual sampling phase rather than in searching and equilibration); (ii) robust (the method is only weakly sensitive to the details of the build‐up protocol, even for relatively short build‐up times); (iii) versatile (a LEUS biasing potential database could easily be preoptimized for small molecules and assembled on a fragment basis for larger ones). © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

GENESIS 1.1: A hybrid‐parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms

GENeralized‐Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all‐atom and coarse‐grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time‐step integration, and hybrid (CPU + GPU) computing. The string method and replica‐exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free‐energy pathway and obtaining free‐energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research. © 2017 Wiley Periodicals, Inc.     

Scalable free energy calculation of proteins via multiscale essential sampling

A multiscale simulation method, "multiscale essential sampling (MSES)," is proposed for calculating free energy surface of proteins in a sizable dimensional space with good scalability. In MSES, the configurational sampling of a full-dimensional model is enhanced by coupling with the accelerated dynamics of the essential degrees of freedom. Applying the Hamiltonian exchange method to MSES can remove the biasing potential from the coupling term, deriving the free energy surface of the essential degrees of freedom. The form of the coupling term ensures good scalability in the Hamiltonian exchange. As a test application, the free energy surface of the folding process of a miniprotein, chignolin, was calculated in the continuum solvent model. Results agreed with the free energy surface derived from the multicanonical simulation. Significantly improved scalability with the MSES method was clearly shown in the free energy calculation of chignolin in explicit solvent, which was achieved without increasing the number of replicas in the Hamiltonian exchange.     

Metadynamics in the conformational space nonlinearly dimensionally reduced by Isomap

Atomic motions in molecules are not linear. This infers that nonlinear dimensionality reduction methods can outperform linear ones in analysis of collective atomic motions. In addition, nonlinear collective motions can be used as potentially efficient guides for biased simulation techniques. Here we present a simulation with a bias potential acting in the directions of collective motions determined by a nonlinear dimensionality reduction method. Ad hoc generated conformations of trans,trans-1,2,4-trifluorocyclooctane were analyzed by Isomap method to map these 72-dimensional coordinates to three dimensions, as described by Brown and co-workers [J. Chem. Phys. 129, 064118 (2008)]. Metadynamics employing the three-dimensional embeddings as collective variables was applied to explore all relevant conformations of the studied system and to calculate its conformational free energy surface. The method sampled all relevant conformations (boat, boat-chair, and crown) and corresponding transition structures inaccessible by an unbiased simulation. This scheme allows to use essentially any parameter of the system as a collective variable in biased simulations. Moreover, the scheme we used for mapping out-of-sample conformations from the 72D to 3D space can be used as a general purpose mapping for dimensionality reduction, beyond the context of molecular modeling.     

On the choice of a reference state for one‐step perturbation calculations between polar and nonpolar molecules in a polar environment

One‐step perturbation is an efficient method to estimate free energy differences in molecular dynamics (MD) simulations, but its accuracy depends critically on the choice of an appropriate, possibly unphysical, reference state that optimizes the sampling of the physical end states. In particular, the perturbation from a polar moiety to a nonpolar one and vice versa in a polar environment such as water poses a challenge which is of importance when estimating free energy differences that involve entropy changes and the hydrophobic effect. In this work, we systematically study the performance of the one‐step perturbation method in the calculation of the free enthalpy difference between a polar water solute and a nonpolar “water” solute molecule solvated in a box of 999 polar water molecules. Both these polar and nonpolar physical reference states fail to predict the free enthalpy difference as obtained by thermodynamic integration, but the result is worse using the nonpolar physical reference state, because both a properly sized cavity and a favorable orientation of the polar solute in a polar environment are rarely, if ever, sampled in a simulation of the nonpolar solute in such an environment. Use of nonphysical soft‐core reference states helps to sample properly sized cavities, and post‐MD simulation rotational and translational sampling of the solute to be perturbed leads to much improved free enthalpy estimates from one‐step perturbation. © 2012 Wiley Periodicals, Inc.     

Replica state exchange metadynamics for improving the convergence of free energy estimates

Metadynamics (MTD) is a powerful enhanced sampling method for systems with rugged energy landscapes. It constructs a bias potential in a predefined collective variable (CV) space to overcome barriers between metastable states. In bias‐exchange MTD (BE‐MTD), multiple replicas approximate the CV space by exchanging bias potentials (replica conditions) with the Metropolis–Hastings (MH) algorithm. We demonstrate that the replica‐exchange rates and the convergence of free energy estimates of BE‐MTD are improved by introducing the infinite swapping (IS) or the Suwa‐Todo (ST) algorithms. Conceptually, IS and ST perform transitions in a replica state space rather than exchanges in a replica condition space. To emphasize this, the proposed scheme is called the replica state exchange MTD (RSE‐MTD). Benchmarks were performed with alanine polypeptides in vacuum and water. For the systems tested in this work, there is no significant performance difference between IS and ST. © 2015 Wiley Periodicals, Inc.     

Free‐energy differences between states with different conformational ensembles

Multiple conformations separated by high‐energy barriers represent a challenging problem in free‐energy calculations due to the difficulties in achieving adequate sampling. We present an application of thermodynamic integration (TI) in conjunction with the local elevation umbrella sampling (LE/US) method to improve convergence in alchemical free‐energy calculations. TI‐LE/US was applied to the guanosine triphosphate (GTP) to 8‐Br‐GTP perturbation, molecules that present high‐energy barriers between the anti and syn states and that have inverted preferences for those states. The convergence and reliability of TI‐LE/US was assessed by comparing with previous results using the enhanced‐sampling one‐step perturbation (OSP) method. A linear interpolation of the end‐state biasing potentials was sufficient to dramatically improve sampling along the chosen reaction coordinate. Conformational free‐energy differences were also computed for the syn and anti states and compared to experimental and theoretical results. Additionally, a coupled OSP with LE/US was carried out, allowing the calculation of conformational and alchemical free energies of GTP and 8‐substituted GTP analogs. © 2013 Wiley Periodicals, Inc.     

Adaptively biased molecular dynamics for free energy calculations

We present an adaptively biased molecular dynamics (ABMD) method for the computation of the free energy surface of a reaction coordinate using nonequilibrium dynamics. The ABMD method belongs to the general category of umbrella sampling methods with an evolving biasing potential and is inspired by the metadynamics method. The ABMD method has several useful features, including a small number of control parameters and an O(t) numerical cost with molecular dynamics time t. The ABMD method naturally allows for extensions based on multiple walkers and replica exchange, where different replicas can have different temperatures and/or collective variables. This is beneficial not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of useful information that may be obtained from a given simulation. The workings of the ABMD method are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated environment.     

Determination of free energy profiles by repository based adaptive umbrella sampling: bridging nonequilibrium and quasiequilibrium simulations

We propose a new adaptive sampling approach to determine free energy profiles with molecular dynamics simulations, which is called as "repository based adaptive umbrella sampling" (RBAUS). Its main idea is that a sampling repository is continuously updated based on the latest simulation data, and the accumulated knowledge and sampling history are then employed to determine whether and how to update the biasing umbrella potential for subsequent simulations. In comparison with other adaptive methods, a unique and attractive feature of the RBAUS approach is that the frequency for updating the biasing potential depends on the sampling history and is adaptively determined on the fly, which makes it possible to smoothly bridge nonequilibrium and quasiequilibrium simulations. The RBAUS method is first tested by simulations on two simple systems: a double well model system with a variety of barriers and the dissociation of a NaCl molecule in water. Its efficiency and applicability are further illustrated in ab initio quantum mechanics/molecular mechanics molecular dynamics simulations of a methyl-transfer reaction in aqueous solution.     

Essential energy space random walks to accelerate molecular dynamics simulations: convergence improvements via an adaptive-length self-healing strategy

Recently, accelerated molecular dynamics (AMD) technique was generalized to realize essential energy space random walks so that further sampling enhancement and effective localized enhanced sampling could be achieved. This method is especially meaningful when essential coordinates of the target events are not priori known; moreover, the energy space metadynamics method was also introduced so that biasing free energy functions can be robustly generated. Despite the promising features of this method, due to the nonequilibrium nature of the metadynamics recursion, it is challenging to rigorously use the data obtained at the recursion stage to perform equilibrium analysis, such as free energy surface mapping; therefore, a large amount of data ought to be wasted. To resolve such problem so as to further improve simulation convergence, as promised in our original paper, we are reporting an alternate approach: the adaptive-length self-healing (ALSH) strategy for AMD simulations; this development is based on a recent self-healing umbrella sampling method. Here, the unit simulation length for each self-healing recursion is increasingly updated based on the Wang-Landau flattening judgment. When the unit simulation length for each update is long enough, all the following unit simulations naturally run into the equilibrium regime. Thereafter, these unit simulations can serve for the dual purposes of recursion and equilibrium analysis. As demonstrated in our model studies, by applying ALSH, both fast recursion and short nonequilibrium data waste can be compromised. As a result, combining all the data obtained from all the unit simulations that are in the equilibrium regime via the weighted histogram analysis method, efficient convergence can be robustly ensured, especially for the purpose of free energy surface mapping.     

Mapping potential energy surfaces

A recently proposed dynamical method [A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002)] allows us to globally sample the free energy surface. This approach uses a coarse-grained non-Markovian dynamics to bias microscopic atomic trajectories. After a sufficiently long simulation time, the global free energy surface can be reconstructed from the non-Markovian dynamics. Here we apply this scheme to study the T=0 free energy surface, i.e., the potential energy surface in coarse-grained space. We show that the accuracy of the reconstructed potential energy surface can be dramatically improved by a simple postprocessing procedure with only minor computational overhead. We illustrate this approach by conducting conformational analysis on a small organic molecule, demonstrating its superiority over traditional unbiased approaches in sampling potential energy surfaces in coarse-grained space.     

A bias-exchange approach to protein folding

By suitably extending a recent approach [Bussi, G.; et al. J. Am. Chem. Soc. 2006, 128, 13435] we introduce a powerful methodology that allows the parallel reconstruction of the free energy of a system in a virtually unlimited number of variables. Multiple metadynamics simulations of the same system at the same temperature are performed, biasing each replica with a time-dependent potential constructed in a different set of collective variables. Exchanges between the bias potentials in the different variables are periodically allowed according to a replica exchange scheme. Due to the efficaciously multidimensional nature of the bias the method allows exploring complex free energy landscapes with high efficiency. The usefulness of the method is demonstrated by performing an atomistic simulation in explicit solvent of the folding of a Triptophane cage miniprotein. It is shown that the folding free energy landscape can be fully characterized starting from an extended conformation with use of only 40 ns of simulation on 8 replicas.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.     

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.     

Free Energy Calculations using a Swarm‐Enhanced Sampling Molecular Dynamics Approach

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm‐enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT‐TI) method. Free energy changes for transitions computed by using IT‐TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm‐enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.     

The reweighted path ensemble

We introduce a reweighting scheme for the path ensembles in the transition interface sampling framework. The reweighting allows for the analysis of free energy landscapes and committor projections in any collective variable space. We illustrate the reweighting scheme on a two dimensional potential with a nonlinear reaction coordinate and on a more realistic simulation of the Trp-cage folding process. We suggest that the reweighted path ensemble can be used to optimize possible nonlinear reaction coordinates.