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 non-equilibrium 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.


INTRODUCTION
When investigating the equilibrium properties of a complex polyatomic system, it is customary to identify a suitable reaction coordinate σ(r 1 , . . . , r N ) : R 3N → Q that maps atomic positions r 1 , . . . , r N onto the points of some manifold Q, and then to study its equilibrium probability density: p(ξ) = δ [ξ − σ(r 1 , . . . , r N )] , ξ ∈ Q (angular brackets denote an ensemble average). The density p(ξ) provides information about the relative stability of states corresponding to different values of ξ along with useful insights into the transitional kinetics between various stable states. In practice, the Landau free energy 1 f (ξ) = −k B T ln p(ξ), is typically preferred over p(ξ), because it tends to be more intuitive. Either p(ξ) or f (ξ) is said to provide a coarse-grained description of the system -in terms of ξ alone -with the rest of the degrees of freedom of the original system integrated out. Quite naturally, the reaction coordinate (often also referred to as collective variable or order parameter ) is typically chosen to represent the slowest degrees of freedom of the original system, although this is not formally required.
In the past few years, several methods targeting the computation of f (ξ) using non-equilibrium dynamics have become popular. First methods that introduced a time evolving potential to bias the original potential energy were the the Local Elevation Method (LEM) 2 , by Huber, Torda and van Gunsteren in the MD context and the Wang-Landau approach in MC one 3 . More recent approaches include the adaptive-force bias method 4 , and the non-equilibrium metadynamics 5,6 method. These methods all estimate the free energy of the reaction coordinate from an "evolving" ensemble of realizations 7,8 , and use that estimate to bias the system dynamics, so as to flatten the effective free energy surface. Collectively, they can all be considered as umbrella sampling methods, with an evolving potential. In the long time limit, the biasing force is expected to compensate for the free energy gradient, so that the biasing potential eventually reproduces the free energy surface.
In this work, we present an Adaptively Biased Molecular Dynamics (ABMD) method whose implementation is particularly efficient and suited for free energy calculations. The method has an O(t) scaling with molecular dynamics time t and is characterized by only two control parameters. In addition, the method allows for extensions based on multiple walkers and replica exchange for both temperature and/or the collective variables. The ABMD method has been implemented in the AMBER software package 9 , and is to be distributed freely.
Before discussing ABMD, it is helpful to review the salient features of the metadynamics (MTD) method. Essentially, the MTD method is built upon the LEM method by exploiting Car-Parrinello dynamics: the phase space of the system is extended to include additional dynamical degrees of freedom harmonically coupled to the collective variable. These additional degrees of freedom are assumed to have masses associated with them, and evolve in time according to Newton's laws. The masses are supposed to be large enough, so that the dynamics of these extra-variables is driven by the free energy gradient. Their trajectory is then used to construct a historydependent biasing potential by means of placing many small Gaussians along the trajectory. When combined with Car-Parrinello ab-initio dynamics, MTD has been successfully used to explore complex reaction pathways involving several energy barriers 10,11,12,13,14,15,16,17,18 .
While MTD continues to be used successfully, there are several known limitations associated with the initial implementation of the method, which provided the motivation for the development of the ABMD method. First, in order to calculate reliable free energies with a controllable accuracy, long runs are needed, especially for the "corrective" follow-up at equilibrium 19 . This is especially true for biomolecular systems, which typically are char-Typeset by REVT E X acterized by many degrees of freedom and non-negligible entropy contributions to the free energies. Long runs, however, may be precluded by the MTD method, because of its unfavorable scaling with molecular dynamics MD time t. While one can readily speed up the original MTD method using such tricks as truncated Gaussians and kd-trees 19 , the bottleneck there is the explicit calculation of the history-dependent potential. Since at every MD step, Gaussians from all previous time steps need to be added, the number of Gaussians grows linearly with t. The numerical cost of MTD therefore grows as O(t 2 ) which, in some cases, may prove itself to be prohibitively expensive, especially when long runs are needed. Another undesirable feature is that the MTD method (at least in its original implementation) is characterized by a relatively large number of parameters (e.g., the masses and spring constants associated with the collective variable, the characteristics of the Gaussians to be added, multiple control parameters, etc), all of which influence the dynamics in an entangled and non-transparent way. A successful MTD run often requires a careful balancing of these parameters, which is especially nontrivial for multidimensional collective variables. More recent implementations of MTD 20 have reduced the number of parameters. As will be discussed, the ABMD method is characterized by only two control parameters and scales as O(t) with simulation time.

II.
THE ADAPTIVELY BIASED MOLECULAR DYNAMICS METHOD.
The ABMD method is formulated in terms of the following set of equations: where the first set represents Newton's equations that govern ordinary MD (temperature and pressure regulation terms are not shown) augmented with the additional force coming from the time-dependent biasing potential U (t|ξ) (with U (t = 0|ξ) = 0), whose time evolution is given by the second equation. In the following, we refer to τ F as flooding timescale and to G(ξ) as kernel (in analogy to the kernel density estimator widely used in statistics 21 ). The kernel is supposed to be positive definite (G(ξ) > 0) and symmetric (G(−ξ) = G(ξ)). It can be perceived as a smoothed Dirac delta function. For large enough τ F and small enough width of the kernel, the biasing potential U (t|ξ) converges towards −f (ξ) as t → ∞ 7,8 . Our numerical implementation of the ABMD method involves the following. We stick with Q = Q 1 × · · · × Q D where Q k is either [a, b] ∈ R 1 or a one-dimensional torus, and use cubic B-splines (or products of thereof for D > 1) to discretize U (t|ξ) in Q: We use the biweight kernel 21 for G(ξ): and an Euler-like discretization scheme for the time evolution of the biasing potential: where σ = σ(r 1 , . . . , r N ) is at time t. Note that this time discretization may be readily improved. This, however, is not really important here, since the goal is not to recover the solution of the ABMD equations per se, but rather to flatten U (t|ξ) + f (ξ) in the t → ∞ limit. Note also, that the numerical cost of evaluation of the time-dependent potential is constant over time, and so ABMD scales trivially as O(t), which is computationally quite favorable. The storage requirements of the ABMD are also quite reasonable, especially if sparse arrays are used for U m . In addition, it is characterized by only two control parameters: the flooding timescale τ F and the kernel width 4∆ξ. ABMD admits two important extensions. The first is identical in spirit to the multiple walkers metadynamics 7,22 . It amounts to carrying out several different MD simulations biased by the same U (t|ξ), which evolves via: where α labels different MD trajectories. A second extension is to gather several different MD trajectories, each bearing its own biasing potential and, if desired, its own distinct collective variable, into a generalized ensemble for "replica exchange" with modified "exchange" rules 23,24,25 . Both extensions are advantageous and lead to a more uniform flattening of U (t|ξ) + f (ξ) in Q. This enhanced convergence to f (ξ) is due to the improved sampling of the "evolving" canonical distribution.
We have implemented the ABMD method in the AMBER package 9 , with support for both replica exchange and multiple-walkers. In pure "parallel tempering" replica exchange (same collective variable in all replicas), N r replicas are simulated at different temperatures T n , n = 1, . . . , N r . Each replica has its own biasing potential U n (t|ξ), n = 1, . . . , N r , that evolves according to its dynamical equation. Exchanges between neighboring replicas are attempted at a prescribed rate, with an exchange probability given by 23 : where E p denotes the atomic potential energy. The biasing potentials are temperature-bound and converge in the t → ∞ limit to the free energies at their respective temperatures.
We have also implemented a more general replica exchange scheme, where different replicas can have different collective variables and/or temperatures, and can experience either an evolving or a static biasing potential (the latter obviously includes the case of U n (t|ξ) = 0). Exchanges between random pairs of replicas are then tried at a prescribed rate. This method is simply a generalization 25 of the "Hamiltonian replica exchange" method described in Ref. 23, and reduces to it when all biasing potentials are static. The big advantage here is that, by using replicas with different collective variables, it is possible to obtain several one-or two-dimensional projections of the free energy surface for the corresponding variables. This is very useful because it not only increases the amount of information that can be gathered from a given simulation, but it also allows for previously obtained information for a collective variable to be used to compute the free energy associated with different variables. For instance, suppose that in the course of a simulation it becomes apparent that one wishes to address additional questions involving different collective variables. Instead of starting from "scratch", one can re-use the already obtained biasing potentials and thereby greatly accelerate the free energy calculation for the new variables. It is also worth noting that for replicas running at the same temperature, the exchange probability does not depend on the atomic potential energies (Eqs. (2.1)-(2.2) above). This implies that the number of replicas needed to maintain acceptable exchange rates can be made independent of the solvent degrees of freedom, provided that one is interested in the properties of the solute only (so that the collective variables do not depend explicitly on the atomic coordinates of the solvent) and that the structure is adequately solvated. This can be exploited to sample a solute with a minimum amount of solvent, and to accelerate the averaging over the solvent degrees of freedom. These last two applications of the general replica exchange method are illustrated in the next section.

III. CASE STUDY: A SHORT PEPTIDE.
To illustrate the method, we have simulated the hydrophobic Ace-GGPGGG-Nme peptide (sketched in Fig.1) in the gas phase and in a solvated environment, using cyclohexane as the explicit solvent. The free energy of this peptide at T = 300 K in gas phase has previously been investigated with the MTD method 19 , and is characterized by a "double-well" structure (see Fig.8), with the wells corresponding to the peptide in a "globular" (left minimum in the Fig.8) and a β-hairpin (right minimum in the  3: RMS error of the free energy over 3.3Å < Rg < 6.3Å for ABMD simulations at T = 300K with τF = 180 ps (1), 360 ps (2), 720 ps (3) and 1440 ps (4).
time step (∆t) is 1 fs for the parallel tempering simulations, and 2 fs otherwise.
The radius of gyration of the heavy atoms was chosen to be the collective variable: Here, R Σ = a (m a /m Σ ) r a is the center of mass, with m Σ = a m a , and the sum runs over all atoms except hydrogen. The initial configuration is the fully unfolded peptide. A reference free energy profile, whose error 19 in the region of interest is less than 0.15 kcal/mol, was computed for benchmarking purposes (see Appendix A for details). As a measure of the RMS free energy error, the following construction was used: accounts for the arbitrary additive constants in the free energies f 1,2 (ξ). Here a = 3.3Å and b = 6.3Å, which correspond to the physical region of interest (see Fig.8). Figure 2 presents the time dependence of the RMS free energy error for the ABMD and reference MTD run 19 . Both simulations have exactly the same kernel width 4∆ξ = 0.25Å, and flooding timescale τ F = 90 ps that corresponds to the a posteriori hills acceptance rate reported in Ref. 19. It is evident that the AMBD run is more accurate than the corresponding MTD run. The ABMD method owes its better convergence to the smoother time evolution of the biasing potential and accurate discretization in Q. The amount of memory used by the ABMD simulation to store U m values was only ≈ 200 × 8 bytes (considering double precision) for roughly 10 8 tiny "hills" that were accumulated by the end of the run. The reference MTD simulation with merely 5 × 10 3 Gaussians required roughly 25 times more memory for the biasing potential (with only the positions of the Gaussians stored explicitly). One can expect, that ABMD will be even more economical when it comes to dealing with multidimensional collective variables, provided that sparse arrays are used for U m with only non-zero elements being stored explicitly.
Although an a priori error estimate for this type of non-equilibrium simulation is really not feasible, it is expected that the error should decrease for increased τ F . This point is illustrated in Fig.3, which shows the error for increasing values of τ F .
In order to decrease the simulation time required for accurate free energy estimates even further, the multiplewalker variation of ABMD proves to be useful. For a moderate number of walkers, the speedup is nearly linear, with an additional increase in accuracy coming from the better sampling of the "evolving" canonical distribution (see Fig.4). Parallel tempering improves both the speed and the accuracy even more. To this end, we first ran ABMD with τ F = 90 ps using 2, 4, 6 and 8 replicas at T = 300 K, 331 K, 365 K, 403 K, 445 K, 492 K, 543 K and 600 K (during equilibrium MD runs, the peptide configuration jumps between the two minima on a picosecond timescale at T = 600K). In all cases, the E RMS was found to be ∼ 1 kcal/mol, or less as t → ∞ (data not shown). Again, the improvement in accuracy stems from the better sampling of the "evolving" canonical distribu- tion. Then, we ran 8 replicas with smaller values τ F and were surprised that the accuracy does not degrade, even for τ F = 11.2 ps as shown in Fig.5. Finally, we turn to aspects related to the general replica exchange method and illustrate its potential. As already noted, by using replicas with different collective variables and swapping these at prescribed rates, it is possible to obtain projections of the free energy surface for the corresponding variables. It is also possible to use previously obtained information with respect to one collective variable to compute the free energy associated with a different variable. For example, suppose that instead of the one-dimensional free energy profile as as a function of R g already discussed, one realizes that what is actually needed is a two-dimensional profile that includes information with respect to the number of O-H bonds along the backbone. The two-dimensional free energy map is computationally quite expensive, but the calculation can be greatly accelerated with the help of the general replica exchange method. We therefore simulated 8 + 1 = 9 replicas. The eight replicas were run at the previously stated temperatures, with each replica biased by a static (not evolving) biasing potential corresponding to the negated free energy associated with the radius of gyration R g at the corresponding temperature, as shown in Fig.8. The additional ninth replica was run at T = 300 K with ABMD flooding in the two collective variables, i.e., R g and the number of O-H bonds along the backbone as given by where r OH is the distance between a pair of hydrogen . In other words, we "re-use" the previously computed free energies for R g to get the free energy in the two-dimensional space (R g , N OH ): the eight first replicas serve as a "sampling enhancement device" for the ninth replica. The calculation is carried out in two stages: a "coarse" stage (15 ns with τ F = 10 ps and 4∆R g = 0.25Å, 4∆N OH = 0.5) followed by a "fine" stage (50 ns with τ F = 100 ps and 4∆R g = 0.1Å, 4∆N OH = 0.25). In both runs exchanges between four randomly chosen pairs of replicas were attempted every 100 fs. The final free energy map is shown in Fig.6. It is clear that this two-dimensional free energy landscape conveys additional information not contained in the one-dimensional free energy plots already discussed. In particular, it allows for a better characterization of the "globular" states of the Ace-GGPGGG-Nme : specifically, it is apparent from the Fig.6 that there are at least two such states with different values of N OH (both correspond to the left minimum in Fig.8). Of course, one could have re-used the information in the one-dimensional R g profiles to include other collective variables, in addition to N OH .
The general replica exchange ABMD may also be advantageous for explicit solvent simulations, which are often notoriously lengthy. Specifically, if one is interested in the solute and the collective variables do not depend on the solvent degrees of freedom, then the number of replicas required to maintain an adequate exchange rate depends only very weakly on the amount of solvent (which must of course be sufficient as to adequately solvate the structure), provided that all the replicas are simulated at the same temperature. This is because the exchange probability does not explicitly depend on the potential energy difference when the temperature of the replicas is the same. While not every choice of collective variables for different replicas will lead to decent exchange rates, one can nevertheless take advantage of this property and use general replica exchange to enhance the sampling in a solvated environment.
In order to demonstrate the method in this regime, we simulated the Ace-GGPGGG-Nme peptide at T = 300 K solvated by 171 cyclohexane (C 6 H 12 ) molecules (the total number of atoms was 3139) under periodic boundary conditions using the General 28 AMBER Force-Field (GAFF) for the solvent. We used a truncated octahedron cell of fixed size (constant volume) that corresponds to the equilibrium density at T = 300 K (the equilibrium density value was obtained from a 10 ns simulation under constant pressure at T = 300 K). The Particle-Mesh Ewald (PME 29 ) method was used for the electrostatic forces, with a 36 × 36 × 36 FFT grid and an 8Å cutoff for the direct sum (same cutoff was used for van der Waals interactions). First, we ran 10 replicas in the "flooding" mode (i.e., under evolving biasing potentials) using as collective variables the distances r CC between the backbone carbons separated by at least 2 amino-acids (there are 10 such distances for Ace-GGPGGG-Nme ). We ran for 5 ns with τ F = 30 ps and 4∆r CC = 1Å, and then for another 5 ns with τ F = 150 ps and 4∆r CC = 0.5Å, attempting exchanges between five randomly selected pairs every 0.5 ps. As expected, at the start of the simulation, the exchange rate was nearly 100% decreasing later as the biasing potentials were built up. By the end of the simulation, when all possible values of the distances had been covered, the exchange rate was very disparate between different pairs of replicas. However, for every replica there was at least one other replica such that the exchange rate between them was reasonable (i.e., the whole simulation did not degenerate into ten different non-interacting trajectories). We then set τ F = ∞ in these 10 replicas, and added an eleventh replica whose collective variable was chosen as the radius of gyration of the heavy atoms. In other words, as before, we use the ten replicas as a "sampling enhancement device" for the last one. We then ran a two-stage flooding scheme: a 5 ns coarser stage, with τ F = 25 ps and 4∆R g = 0.25Å; followed by a 10 ns finer stage, with τ F = 180 ps and 4∆R g = 0.2Å. As before, the exchange attempts between 5 randomly selected pairs of replicas were performed every 0.5 ps. The "raw" ABMDcomputed free energy associated with R g after that stage is shown in Fig.7. In a next step, we ran 64 biased simulations (each comprising of 11 replicas) for 7 ns (first 2 ns for equilibration followed by 5 ns of "production" runs) starting from different initial configurations. We set τ F = ∞ in all replicas (static biasing potentials) and recorded the values of R g in the eleventh replica every 10 picoseconds. We then used the log-spline algorithm of Stone 30 at. al. to estimate the (biased) log-density of the R g values at equilibrium. This led us to the final shape of the free energy curve shown in Fig.7. Compared to the gas phase, the folded β-turn in the cyclohexane solvated peptide is clearly favored over the globular structure.

IV.
CONCLUSIONS AND OUTLOOK.
In summary, we have presented an ABMD method that computes the free energy surface of a reaction coordinate using non-equilibrium dynamics. The method belongs to the general category of umbrella sampling methods with an evolving potential, and is characterized by only two control parameters (the flooding timescale and the kernel width) and a favorable O(t) scaling with molecular dynamics time t. This scaling can be very important for large-scale classical MD biomolecular simulations when long simulation times are required (see, for example Ref. 31, and references therein).
ABMD has also been extended to support multiple walkers and replica exchange. Both variations improve speed and accuracy of the method due to the better sampling of the "evolving" canonical distribution. The replica exchange ABMD has been generalized to include different temperatures and/or collective variables, that move under either an evolving or a static biasing potential. Aside from enhancing the sampling, this swapping of replicas has several important practical advantages. Most importantly, it enables one to obtain projections of the free energy surface for any number of collective variables one might wish to investigate. In addition, one can re-use previously obtained results in order to enhance the sampling of new collective variables. It is also possible to exploit the fact that exchange rates at the same temperature are independent of the potential energy to enhance sampling of a solute in a minimum amount of solvent (for collective variables independent of solvent atom coordinates). We have implemented the ABMD method in the AMBER package 9 , and plan to distribute it freely. Here, we have demonstrated the workings of the ABMD method with a study of the folding of the Ace-GGPGGG-Nme peptide, The application of ABMD to more complicated biomolecular systems is reserved for future publications.