Driven Metadynamics: Reconstructing Equilibrium Free Energies from Driven Adaptive-Bias Simulations.

We present a novel free-energy calculation method that constructively integrates two distinct classes of nonequilibrium sampling techniques, namely, driven (e.g., steered molecular dynamics) and adaptive-bias (e.g., metadynamics) methods. By employing nonequilibrium work relations, we design a biasing protocol with an explicitly time- and history-dependent bias that uses on-the-fly work measurements to gradually flatten the free-energy surface. The asymptotic convergence of the method is discussed, and several relations are derived for free-energy reconstruction and error estimation. Isomerization reaction of an atomistic polyproline peptide model is used to numerically illustrate the superior efficiency and faster convergence of the method compared with its adaptive-bias and driven components in isolation.

Computational investigation of complex systems such as biological macromolecules is often hampered by the need to sample rare events, that is, crossing high free-energy barriers, which cannot be accomplished by standard sampling techniques such as conventional room-temperature molecular dynamics (MD) simulations. To enhance the sampling of such rare barrier-crossing events, one may use a biasing potential. Using a time-independent bias (e.g., umbrella sampling technique[1]) has the advantage of keeping the system in equilibrium; however, designing a “practical” biasing potential may require a priori knowledge of the free-energy landscape, particularly when this landscape is rugged. On the other hand, simulations generated by time-dependent biasing protocols require a nonequilibrium treatment. In the context of MD simulations using (i) a history-dependent “adaptive-bias” force/potential and (ii) an external (explicitly time-dependent) “driving” force/potential are two well-known[2−6] nonequilibrium schemes to enhance the sampling and estimate the free energies.

In adaptive-bias methods, for example, local elevation,[2] coarse MD,[7] adaptively biased MD,[8] adaptive accelerated MD,[9] metadynamics,[3,4] and λ-metadynamics,[10] the simulation history is used to enhance the sampling by discouraging the system to return to the already visited regions of the phase space. Potential is adaptively biased until converged and used to reconstruct the free energy landscape. A nonequilibrium driven scheme (e.g., steered MD[11]) can be used to induce a transition by steering the system along a reaction coordinate that otherwise may not be sampled in an unbiased simulation. By measuring the work along these nonequilibrium trajectories, one may estimate the free-energy differences or reconstruct the free-energy profiles using nonequilibrium work relations.[5,6,12] These relations can be also used to estimate relative reaction rates[13] and generally any equilibrium path-ensemble averages.[14]

Adaptive-bias and driven MD schemes are powerful methods whose applications have been extended to biomolecular simulations;[15−19] however, both methods have practical limitations. Steered MD is often associated with a slow convergence if used for free-energy calculations (because the dissipative work is often large and many iterations may be needed to capture a small one, that is, the dominant term in work-based free energy estimators[5]) but it can be used to explore the transition paths, at least qualitatively, an advantage over metadynamics, in which the system starting from one end of the configuration space (the reactant) may take a long time to visit the other end (the product). One cannot estimate the free-energy difference of the two end states until both states have been sufficiently sampled.

By incorporating an explicitly time-dependent scheme into a history-dependent scheme, we introduce a novel driven adaptive-bias scheme, termed driven metadynamics (d-MetaD), that takes advantage of both its driven and adaptive-bias components and is advantageous over both components in isolation. d-MetaD has an advantage over conventional metadynamics in that it ensures the exploration of the transition pathway (from one end to the other) in the early stages of the simulation and gradually improves the estimate of the free energies almost uniformly along the reaction coordinate. d-MetaD also has an advantage over the conventional steered MD in that the effective free-energy surface gradually becomes smooth and flat such that the system can move along the reaction coordinate with progressively less amount of work.

We use nonequilibrium work relations not only in analyzing the trajectories but also more importantly in constructing the biasing protocol that is per se a novel use of nonequilibrium work relations. The approach presented here expands the scope of nonequilibrium work relations into a new realm, namely, designing work-based biasing protocols. We introduce a general framework to combine an adaptive-bias scheme with a nonequilibrium driven scheme by reweighting the kernel memory (used in most adaptive-bias schemes) using on-the-fly nonequilibrium work measurements. Although, for the sake of clarity, we discuss a particular algorithm here, this approach can be easily generalized to combine any methods from the adaptive-bias class with any from the nonequilibrium driven class, resulting in “hybrid” protocols that could be more efficient than their adaptive-bias and driven components in isolation.

Consider a system described by coordinates r and momenta p and governed by the Hamiltonian H(r,p) with a canonical equilibrium distribution at temperature 1/β. Suppose that x(r) is a holonomic coordinate (e.g., representing a slow mode) and F(x) is the free energy associated with x; that is, p(x) = e–β describes the distribution of x in equilibrium:

In the adaptive-bias scheme, one adds a history-dependent biasing potential to the system such as:in which K(Δx) is a kernel function (e.g., (1/((2π)1/2σ))exp(−(1/2)(Δx/σ)2) and σ is the kernel width), U0(x) is an arbitrary function, and ω(x, t) is an energy rate that is ω0e–β in the well-tempered metadynamics[20] (1/β′ is an arbitrary temperature) that reduces to a constant ω0 in the β′ → 0 limit (i.e., conventional metadynamics[3]). Kernel K(Δx) approximates a smooth δ function, and ⟨K(x – x)⟩acan be thought of as an approximation of the probability density p(x, t) (⟨.⟩a is the average over an ensemble of adaptive-bias trajectories). It has been shown[21,20] that ⟨Ua(x,t → ∞)⟩a ≈ Uas(x) + u(t) (u(t) is an additive constant at each t) in which Uas(x) equals −F(x) and −(1 + (β′/β))−1F(x) in the conventional[21] and well-tempered[20] metadynamics, respectively. The steady-state distribution of x can be described by ps(x) ∝ e–β(. U0(x) could be a flat function or an initial guess for Uas(x) to speed up the convergence process, but Uas(x) does not depend on it. Here we focus only on these variations, although the generalization of our method to other variations of the adaptive-bias scheme is straightforward.

In the nonequilibrium driven scheme, one adds a driving potential Ud(x, t) to the system that is often harmonic, Ud(x, t) = (1/2)k(x – X(t))2, in which X(t) is a protocol controlling the target (e.g., X(t) = x0 + (x1 – x0)(t/T); x0 and x1 may represent the reactant and product or the two ends of the x space to be explored). Although the system might remain far from equilibrium, one can reconstruct the unbiased distribution[6] according to:in which Δw = w – Ud(x,t), w = ∫0 dt′(∂/∂t′)Ud(x,t′) is the total work up to time t, and ⟨.⟩d denotes an average over the ensemble of driven trajectories (driven ensemble). The initial configurations are prepared in an equilibrium state perturbed by the biasing potential Ud(x,0). The e–βΔ term can be thought of as a weighting factor to connect the driven ensemble to the equilibrium ensemble.[22]

To combine the two schemes described above, here we introduce a driven adaptive-bias scheme that adds an adaptive (Ua(x, t)) and a driving (Ud(x, t)) potential to the Hamiltonian. We use an iterative approach in which an independent simulation is performed from time t = 0 to T in the nth iteration (n = 1, 2, ...), biased by the potential Ud(x, t) + Ua(x, t) in which Ud(x, t) = (k/2)(x – X(t))2 for all n and:in which Δw = w – Ud(x, t) similar to relation 3 and U(x) = Ua(x, T) for n > 0 and U0(x) is an initial guess. ω(x, t) is the same as its nondriven (conventional or well-tempered) counterpart. Note that the “weighted” average ⟨K(x – x)e–βΔ⟩ approximates ⟨δ(x – x)e–βΔ⟩da (⟨.⟩da denotes an average over the driven adaptive-bias ensemble). We also note that for practical reasons the e–βΔ factor in relation 4 can be also reweighted using “weight functions” similar to those used for the free-energy estimation from conventional pulling experiments. (See ref (23).) To simplify the discussions we leave the e–βΔ factor as is (i.e., using the “constant weight” protocol[23]). However, the “pulling potential” or “occupancy” weights are generally more appropriate for practical/numerical reasons.[23]

In general, the initial configuration for iteration n needs to be equilibrated under the biasing potential Ua(x,0) + Ud(x,0). However, preparing the initial configurations can be simplified if we assume that the driven term is considerably greater than the adaptive term, (guaranteed in the stiff-spring limit[24]) such that the equilibrium distribution for the perturbed system is more or less independent of the adaptive term. Because the driven potential is the same for all n, the initial configurations for all iterations can be also prepared simply by using the Ud(x,0) bias.

If we iterate the protocol described above until the adaptive term of the biasing potential converges to Uas(x), the system can be considered as a nonequilibrium driven system governed by the Hamiltonian H(r,p) + Uas(x) + Ud(x, t). One can show that Uas(x) of a driven metadynamics protocol is the same as that of its nondriven counterpart (see Supporting Information), for example, Uas(x) = −(1 + (β′/β))−1F(x) for the well-tempered variation and Uas(x) = −F(x) for the β′ → 0 limit.

Although Uas(x) can be used directly to estimate the free energy of the system, one might alternatively measure work along the trajectories once the Ua(x, t) has nearly converged (i.e., Ua(x, t) ≈ Uas(x) + u(t)). If the driving process at time t = 0 starts from the state X(0), governed by the Hamiltonian H(r,p) + Ua(x, t) + Ud(x,0), then the distribution of states associated with x in the steady state, governed by H(r,p) + Ua(x, t), may be represented by ps(x) ∝ exp(−βFs(x)), and Fs(x) = F(x) + Uas(x) may be estimated from the nonequilibrium work relation:

In principle, this estimator will result in an accurate description of free energy along x by using the work measurements at any given t, assuming enough sampling for all x. In practice, one might combine multiple time slices (by storing x and w at different t) to estimate the free energy, for instance, by employing the Hummer–Szabo estimator method:[6]that can be solved self-consistently along with f(t):Ideally, in the case of conventional variation of metadynamics, Fs(x) is flat because Uas(x) = −F(x). In practice, Fs(x) estimated from the relation above is not flat and can be used to correct the free energy estimated from biasing potential FU(x) = −Uas(x) via FW = FU(x) + Fs(x); FW is the corrected free energy, and ε(x) = |FW(x) – FU(x)| = |Fs(x)| is an error estimate associated with it. This can be thought of as the driven counterpart of “umbrella corrections”[8] in conventional adaptive-bias simulations. In the case of well-tempered metadynamics, Fs(x) = −(β′/β)F(x), so the above work-based estimate of F(x) can be used as an alternative method for estimating F(x) as −(β/β′)Fs(x). However, choosing a β′ considerably smaller than β may result in a larger error than direct use of Uas(x). In this case, Fs(x) can be used to correct FU(x) = −(1 + (β/β′))Uas(x) and estimate the error similar to the case of β′ = 0 discussed above. The corrected free energy would be FW(x) = −Uas(x) + Fs(x), associated with an error ε(x) = |FW(x) – FU(x)| = |Fs(x) + (β′/β)Uas(x)|. One can estimate Fs(x) and Uas(x) from work measurements and converged potential, respectively.

We note that one may derive other work-based free-energy estimators for driven adaptive-bias trajectories similar to those used in nonadaptive variations of nonequilibrium driven systems. More importantly, a bidirectional driving protocol is known to result in a faster convergence than a unidirectional one by combining the forward and reverse driven ensembles. d-MetaD as introduced here can be used in a bidirectional scheme such that the forward direction XF(t) varies from XF(0) to XF(T) and in the reverse direction XR(t) = XF(T – t). One can iteratively run the forward and reverse simulations in an alternating manner with an initial equilibration for each simulation. Note that the two protocols (forward and reverse) are associated with the same Uas(x). Deriving bidirectional free-energy estimators[12] is also straightforward.

d-MetaD in ω0 → 0 and k → 0 limits describes the nonadaptive driven MD and nondriven metadynamics, respectively. Interestingly, because of its iterative nature, at the end of any iteration n, d-MetaD can be easily turned into metadynamics or steered MD by removing its driven or adaptive component, respectively. The biasing potential U(x) will be the initial guess for the metadynamics simulations and a time-independent bias for the steered MD simulations. By removing both driven and adaptive components and using U(x) as a time-independent bias, d-MetaD can be turned into a conventional biased equilibrium simulation.

Finally, our method can be generalized to reconstruct multidimensional free-energy surfaces. While this generalization is straightforward in adaptive-bias simulations, the driving protocol is easier to implement using a 1-D reaction coordinate, x. One can use the driving potential Ud(x, t) and construct a multidimensional adaptive potential Ua(x,y,t) described by:in which y represents collective variables other than x and K′ is a multidimensional kernel function.

To test the d-MetaD method and compare its performance to that of metadynamics and steered MD, we carried out several sets of MD simulations, investigating the cis–trans isomerization reaction in polyproline peptides. Proline is unique among the common amino acids in that its side chain is cyclized onto the backbone nitrogen atom, restricting its conformational space to a narrow region in the Ramachandran map.[25] However, of the 20 common amino acids, only proline is “comfortable” in both cis- and trans-isomer conformations. Polyproline is known to form helical structures with two well-characterized conformations: (i) a right-handed polyproline type-I helix (PPI) with all residues in the cis-isomer conformation (all-cis conformer) and (ii) a left-handed polyproline type II helix (PPII) with all residues in the trans-isomer conformation (all-trans conformer) (see Figure 1). Traditionally, an oligomeric proline peptide in aqueous solution is seen as a “rigid rod”; however, recent experimental[26] and computational[16] studies suggest a considerable heterogeneity in the conformational space of polyproline with subpopulations of distinct end-to-end distances due to the existence of distinct cis–trans patterns.

Figure 1

Top (ribbon representation) and side (licorice representation) views of the right-handed PPI and the left-handed PPII conformations of a pentameric polyproline peptide. The backbone atoms involved in the definition of ω dihedral angles are shown; the prolyl peptide bonds are highlighted.

From a computational point of view, the characterization of proline-rich peptides is rather difficult because the cis–trans isomerization reaction is an extremely slow process (i.e., tens to hundreds of seconds at room temperature[27]). This reaction has recently been studied using conventional adaptive-bias and driven MD schemes both on pure polyproline[16,28] and on other proline-containing peptides.[29−31] Free-energy landscapes of polyproline peptides in different environments including water, propanol, and hexane, obtained using a variation of conventional metadynamics (i.e., adaptively biased MD[8]) reveal that water and hexane favor PPII while propanol favors PPI.[16,28] Here we investigate the PPII ↔ PPI reaction of polyproline peptides (Ac-(Pro)-NH2) of various length (1 ≤ n ≤ 5) in an implicit water environment (using the Generalized Born model[32]). NAMD 2.8[33] simulation package was used to perform the simulations using the CHARMM27 force field[34] with a 1 fs time step at a constant temperature of 300 K (using a Langevin thermostat with a damping coefficient of 1/ps) with no cutoff for nonbonded interactions.

To study the cis–trans isomerization, we define a collective variable Ω = ∑cos2(ω/2) based on the backbone dihedral angles {ω} that roughly measures the number of cis-prolyl bonds in an n-mer peptide; ω is around 0 (±180) for cis (trans) isomers, and thus cos2(ω/2) falls around 1 (0). Ω ranges between 0 and n (representing all-trans PPII and all-cis PPI conformers, respectively) and is associated with a free energy F(Ω) with n + 1 minima and n maxima centered around integer and half-integer values of Ω, respectively.[16,28] We note that generally designing a 1D collective variable (e.g., Ω) that reflects all of the slow conformational changes of the system (e.g., cis/trans isomerization of all prolyl bonds) is a key factor in free-energy calculation methods such as metadynamics, steered MD, and inherently d-MetaD.

A bidirectional steered MD (SMD) scheme was used to steer the system along Ω between 0 and n with a constant speed of Ω̇ = 10/ns and a harmonic constant of 100 kcal/mol, iterated 10 000 times in each direction with a total simulation time of 2, 4, 6, 8, and 10 μs for monomeric, dimeric, trimeric, tetrameric, and pentameric polyproline peptides, respectively. A metadynamics MD (MetaD) scheme described by a Gaussian kernel of width σ = 0.2 and a well-tempered rate of (ω0 = 0.1 (kcal/mol) ps–1, (β/β′) = 10), and a flat initial bias U0(Ω) = 0 was used to reconstruct F(Ω) using a total simulation time matching that of the SMD runs for each peptide. Finally, a d-MetaD scheme combining the SMD and MetaD methods with the exact same parameters described above was used to reconstruct F(Ω). The e–βΔ factor was simplified to e–β, mimicking the pulling potential weight protocol[23] (i.e., e–βΔ was multiplied by the unnormalized pulling potential weight e–β).

In the SMD and d-MetaD simulations, Ω and work values were stored every 1 ps and used at the end of each bidirectional iteration to reconstruct the free-energy profiles F(Ω) and Fs(Ω) (for SMD and d-MetaD, respectively) using the self-consistent bidirectional implementation[12] of the Hummer–Szabo method.[6] For the MetaD and d-MetaD simulations, the instantaneous adaptive potential Ua(Ω,t) was used as an estimate for Uas(Ω) to reconstruct F(Ω) at time t. The corrected d-MetaD free energy FW(Ω) was constructed using Uas(Ω) and Fs(Ω), as approximated from Ua(Ω,t) and work measurements, respectively.

Figure 2 summarizes our results for the pentameric peptide; SMD, MetaD, and d-MetaD methods give very similar estimates for F(Ω) eventually (Figure 2a), but the d-MetaD simulations converge considerably faster than the others (see Figure 2b). To illustrate the convergence behavior, we plot two physically relevant quantities in Figure 2c,d: the free-energy difference between the PPI and the PPII conformers, F(5) – F(0), and between the transition state and the global minimum, F(2.5) – F(0). Figure 2 generally reveals that d-MetaD constructively combines both MetaD and SMD methods to speed up the convergence time. To quantify the convergence time associated with a method m (SMD, MetadD, or d-MetaD), we define it as the smallest time τc in which |F(Ω) – Fc(Ω)| < ε for all Ω and for all t > τc; F(Ω) is the free energy estimated from the method m at time t, Fc(Ω) is our final free-energy estimate averaged over the results of all three methods, and ε is an error parameter adjusted to ensure the meaningfulness of the analysis statistically. Figure 3 shows τc associated with all three methods for all of the peptides studied using ε = 2 kcal/mol. τc is expected to grow exponentially by the peptide length n (that somewhat represents the complexity of the problem) due to the exponential growth of the number of states; considering the cis-trans isomerization only, we precisely have 2 conformers for an n-mer peptide. Interestingly, the τc growth constant of d-MetaD method is considerably smaller than that of MetaD and SMD such that at least in this particular case and within this range (1 ≤ n ≤ 5) τc grows almost linearly by n.

Figure 2

(a,b) Free-energy profile F(Ω) (offset by F(0)) of a pentameric polyproline peptide obtained from SMD (red), MetaD (blue), and d-MetaD (black) simulations at t = 10 000 and 200 ns, respectively. In panel b, the dashed curve is the d-MetaD results without the work-based corrections (i.e., F(Ω)) and the green curve is the average of the three converged curves in (a) (i.e., F(Ω)). (c,d) Evolution of ΔF = F(5) – F(0) and ΔF = F(2.5) – F(0) by time as estimated from SMD (red), MetaD (blue), and d-MetaD (black). The green regions represent F(Ω) ± ε with ε = 2 kcal/mol.

Figure 3

Convergence time τc associated with SMD (red squares), MetaD (blue circles), and d-MetaD (black triangles) simulations performed on polyproline peptides of length n = 1, ..., 5. The data were fitted using a τc = reα function (linear fitting in the (log τc,n) space). Inset: τc associated with the d-MetaD method fitted using a linear function τc = an.

In summary, we have introduced a nonequilibrium free-energy method with an explicitly time- and history-dependent biasing protocol. d-MetaD is a novel scheme that takes advantage of both driven and adaptive-bias schemes. Unlike the nondriven metadynamics, d-MetaD samples the reaction path almost uniformly from the very beginning, and unlike the nonadaptive driven MD, the amount of work along the reaction path decreases gradually. The history-dependent term of the potential on average behaves similarly to that of conventional adaptive-bias systems and once converged can be used to reconstruct the free-energy profile whose error can be estimated using nonequilibrium work measurements.