Adaptively Biased Molecular Dynamics for Free Energy Calculations

0.35 MB·English

by  Volodymyr Babin

Adaptively Biased Molecular Dynamics for Free Energy Calculations
Volodymyr Babin, Christopher Roland, and Celeste Sagui∗
Center for High Performance Simulations (CHiPS) and Department of Physics,
North Carolina State University, Raleigh, NC 27695
(Dated: January 16, 2008)
WepresentanAdaptivelyBiasedMolecularDynamics(ABMD)methodforthecomputationofthefree
energy surface of a reaction coordinate using non-equilibriumdynamics. The ABMD method belongs
to the general category of umbrella sampling methods with an evolving biasing potential, and is
8 inspired by the metadynamics method. The ABMD method has several useful features, including a
0 small number of control parameters, and an O(t) numerical cost with molecular dynamics time t.
0 The ABMD method naturally allows for extensions based on multiple walkers and replica exchange,
2 wheredifferentreplicascanhavedifferenttemperaturesand/orcollectivevariables. Thisisbeneficial
not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of
n usefulinformationthatmaybeobtainedfromagivensimulation. TheworkingsoftheABMDmethod
a
are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated
J
environment.
4
2

I. INTRODUCTION they can all be considered as umbrella sampling meth-
]
i ods, with an evolving potential. In the long time limit,
c
s When investigating the equilibrium properties of a the biasing force is expected to compensate for the free
- complex polyatomic system, it is customary to identify energy gradient, so that the biasing potential eventually
l
r a suitable reaction coordinate σ(r ,...,r ) : R3N 7→ Q reproduces the free energy surface.
mt that maps atomic positions r1,..1.,rN oNnto the points Inthiswork,wepresentanAdaptivelyBiasedMolecular
of some manifold Q, and then to study its equilibrium Dynamics (ABMD) method whose implementation is par-
.
t probability density: ticularly efficient and suited for free energy calculations.
a
m ThemethodhasanO(t)scalingwithmoleculardynamics
p(ξ)= δ[ξ−σ(r1,...,rN)] , ξ ∈Q time t and is characterized by only two control parame-
-
d (cid:10) (cid:11) ters. Inaddition,themethodallowsforextensionsbased
n (angularbracketsdenoteanensembleaverage). Theden- onmultiplewalkersandreplicaexchangeforbothtemper-
o sityp(ξ)providesinformationabouttherelativestability ature and/or the collective variables. The ABMD method
c ofstatescorrespondingtodifferentvaluesofξ alongwith has been implemented in the AMBER software package9,
[
usefulinsights into the transitionalkinetics betweenvar- and is to be distributed freely.
2 ious stable states. In practice, the Landau free energy1 Before discussing ABMD, it is helpful to review the
v
salient features of the metadynamics (MTD) method. Es-
3 f(ξ)=−k T lnp(ξ),
B sentially, the MTD method is built upon the LEM method
5
4 is typically preferred over p(ξ), because it tends to be by exploiting Car-Parrinello dynamics: the phase space
0 more intuitive. Either p(ξ) or f(ξ) is said to provide of the system is extended to include additional dynami-
8. a coarse-grained description of the system — in terms cal degrees of freedom harmonically coupled to the col-
0 of ξ alone — with the rest of the degrees of freedom of lective variable. These additional degrees of freedom are
7 the original system integrated out. Quite naturally, the assumedtohavemassesassociatedwiththem,andevolve
0 reaction coordinate (often also referred to as collective in time according to Newton’s laws. The masses are
v: variable ororder parameter) istypicallychosento repre- supposed to be large enough, so that the dynamics of
i senttheslowestdegreesoffreedomoftheoriginalsystem, these extra-variables is driven by the free energy gradi-
X
although this is not formally required. ent. Theirtrajectoryisthen usedto constructahistory-
ar In the past few years, several methods targeting the dependent biasing potential by means of placing many
computation of f(ξ) using non-equilibrium dynamics small Gaussians along the trajectory. When combined
have become popular. First methods that introduced with Car-Parrinello ab-initio dynamics, MTD has been
a time evolving potential to bias the original potential successfully used to explore complex reaction pathways
energy were the the Local Elevation Method (LEM)2, by involving several energy barriers10,11,12,13,14,15,16,17,18.
Huber, Torda and van Gunsteren in the MD context and While MTD continues to be used successfully, there are
the Wang-Landau approachin MC one3. More recent ap- several known limitations associated with the initial im-
proaches include the adaptive-force bias method4, and plementationofthe method, whichprovidedthe motiva-
the non-equilibrium metadynamics5,6 method. These tion for the development of the ABMD method. First, in
methods all estimate the free energy of the reaction co- ordertocalculatereliablefreeenergieswithacontrollable
ordinate from an “evolving” ensemble of realizations7,8, accuracy, long runs are needed, especially for the “cor-
andusethatestimatetobiasthe systemdynamics,soas rective” follow-up at equilibrium19. This is especially
to flatten the effective free energy surface. Collectively, true for biomolecular systems, which typically are char-
Typeset by REVTEX
2
acterizedby manydegreesof freedomand non-negligible andusecubicB-splines(orproductsofthereofforD >1)
entropy contributions to the free energies. Long runs, to discretize U(t|ξ) in Q:
however, may be precluded by the MTD method, because
of its unfavorable scaling with molecular dynamics MD U(t|ξ)= U (t)B(ξ/∆ξ−m),
m
time t. While one can readily speed up the original MTD mX∈ZD
method using such tricks as truncated Gaussians and
kd-trees19, the bottleneck there is the explicit calcula-
(2−|ξ|)3/6, 1≤|ξ|<2,
tion of the history-dependent potential. Since at every
MD step, Gaussians from all previous time steps need to B(ξ)= ξ2(|ξ|−2)/2+2/3, 0≤|ξ|<1,
 0, otherwise.
be added, the number of Gaussians grows linearly with
t. The numerical cost of MTD therefore grows as O(t2) 
We use the biweight kernel21 for G(ξ):
which,insomecases,mayproveitselftobeprohibitively
expensive,especiallywhenlongrunsareneeded. Another 48 1− ξ2 4 2, −2≤ξ ≤2
undesirable feature is that the MTD method (at least in G(ξ)= ,
41(cid:26)0(cid:0), (cid:14) (cid:1) otherwise
its original implementation) is characterized by a rela-
tively large number of parameters (e.g., the masses and
and anEuler-like discretizationscheme for the time evo-
spring constants associated with the collective variable, lution of the biasing potential:
the characteristics of the Gaussians to be added, multi-
ple control parameters, etc), all of which influence the kBT
U (t+∆t)=U (t)+∆t G σ/∆ξ−m ,
dynamics in an entangled and non-transparent way. A m m τ
successful MTD run often requires a careful balancing of F (cid:2) (cid:3)
whereσ =σ(r ,...,r )isattimet. Notethatthistime
theseparameters,whichisespeciallynontrivialformulti- 1 N
discretizationmaybereadilyimproved. This,however,is
dimensional collective variables. More recent implemen-
tationsof MTD20 havereducedthenumberofparameters. notreallyimportanthere,sincethegoalisnottorecover
As will be discussed, the ABMD method is characterized the solution of the ABMD equations per se, but rather to
by only two control parameters and scales as O(t) with flatten U(t|ξ)+f(ξ) in the t→∞ limit. Note also, that
the numerical cost of evaluation of the time-dependent
simulation time.
potential is constant over time, and so ABMD scales triv-
ially as O(t), which is computationally quite favorable.
II. THE ADAPTIVELY BIASED MOLECULAR The storage requirements of the ABMD are also quite rea-
DYNAMICS METHOD. sonable, especially if sparse arrays are used for Um. In
addition, it is characterizedby only two control parame-
The ABMDmethod isformulatedintermsofthe follow- ters: thefloodingtimescaleτF andthekernelwidth4∆ξ.
ABMD admits two important extensions. The first is
ing set of equations:
identical in spirit to the multiple walkers metadynam-
d2r ∂ ics7,22. It amounts to carrying out several different MD
ma dt2a =Fa− ∂r U t|σ(r1,...,rN) , simulationsbiasedbythesameU(t|ξ),whichevolvesvia:
a (cid:2) (cid:3)
∂U(t|ξ) k T
= B G ξ−σ(rα,...,rα) ,
∂U(t|ξ) = kBTG ξ−σ(r ,...,r ) , ∂t τF Xα (cid:2) 1 N (cid:3)
1 N
∂t τ
F (cid:2) (cid:3) where α labels different MD trajectories. A second exten-
where the first set represents Newton’s equations that sion is to gather several different MD trajectories, each
govern ordinary MD (temperature and pressure regula- bearing its own biasing potential and, if desired, its
tiontermsarenotshown)augmentedwiththeadditional own distinct collective variable, into a generalized en-
force coming from the time-dependent biasing potential semble for “replica exchange” with modified “exchange”
U(t|ξ) (with U(t = 0|ξ) = 0), whose time evolution is rules23,24,25. Both extensions are advantageous and lead
given by the second equation. In the following, we re- to a more uniform flattening of U(t|ξ)+f(ξ) in Q. This
fer to τ as flooding timescale and to G(ξ) as kernel (in enhancedconvergencetof(ξ)isduetotheimprovedsam-
F
analogy to the kernel density estimator widely used in pling of the “evolving” canonical distribution.
statistics21). The kernel is supposed to be positive defi- We have implemented the ABMD method in the AMBER
nite (G(ξ) > 0) and symmetric (G(−ξ) = G(ξ)). It can package9, with support for both replica exchange and
be perceived as a smoothed Dirac delta function. For multiple-walkers. In pure “parallel tempering” replica
large enough τ and small enough width of the kernel, exchange (same collective variable in all replicas), N
F r
the biasing potential U(t|ξ) converges towards −f(ξ) as replicas are simulated at different temperatures T , n =
n
t→∞7,8. 1,...,N . Each replica has its own biasing potential
r
OurnumericalimplementationoftheABMDmethodin- Un(t|ξ), n= 1,...,N , that evolves according to its dy-
r
volves the following. We stick with Q = Q ×···×Q namicalequation. Exchangesbetweenneighboring repli-
1 D
whereQ iseither [a,b]∈R1 oraone-dimensionaltorus, casareattemptedataprescribedrate,withanexchange
k
3
6
O
NH
5
N
MTD
)
NH ol4
m
O O O al/
O HN kc3
(
S
NH NH M2
R
E
NH 1
ABMD
O O
0
0 50 100 150 200 250 300
FIG.1: TheAce-GGPGGG-Nme peptideinaβ-hairpin confor- time (ns)
mation (sketch).
FIG.2: RMS errorof thefreeenergy over3.3˚A<Rg <6.3˚A
at T = 300K for the ABMD (blue) and reference MTD (red)
probability given by23: simulations.
1, ∆≤0,
w(m|n)= (2.1)
(cid:26)exp(−∆), ∆>0, 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
1 1
∆ = − Em−En (2.2) ontheatomicpotentialenergies(Eqs.(2.1)-(2.2)above).
(cid:18)kBTn kBTm(cid:19)(cid:16) p p(cid:17) Thisimpliesthatthenumberofreplicasneededtomain-
1
+ Um(ξn)−Um(ξm) tainacceptableexchangeratescanbe madeindependent
kBTmh i of the solvent degrees of freedom, provided that one is
1 interested in the properties of the solute only (so that
− Un(ξn)−Un(ξm) ,
k T the collective variables do not depend explicitly on the
B nh i
atomiccoordinatesofthesolvent)andthatthestructure
where E denotes the atomic potential energy. The bi-
p is adequately solvated. This can be exploited to sam-
asing potentials are temperature-bound and converge in
ple a solute with a minimum amount of solvent, and to
the t → ∞ limit to the free energies at their respective
accelerate the averagingover the solvent degrees of free-
temperatures.
dom. These last two applications of the general replica
We have also implemented a more general replica ex-
exchange method are illustrated in the next section.
changescheme,wheredifferentreplicascanhavedifferent
collectivevariablesand/ortemperatures,andcanexperi-
ence either an evolving or a static biasing potential (the
latter obviously includes the case of Un(t|ξ) = 0). Ex- III. CASE STUDY: A SHORT PEPTIDE.
changes between random pairs of replicas are then tried
at a prescribed rate. This method is simply a general- To illustrate the method, we have simulated the hy-
ization25 of the “Hamiltonian replica exchange” method drophobic Ace-GGPGGG-Nme peptide (sketched in Fig.1)
described in Ref.23, and reduces to it when all biasing inthegasphaseandinasolvatedenvironment,usingcy-
potentials are static. The big advantage here is that, clohexane as the explicit solvent. The free energy of this
by using replicas with different collective variables, it is peptide at T = 300K in gas phase has previously been
possibletoobtainseveralone-ortwo-dimensionalprojec- investigatedwiththeMTDmethod19,andischaracterized
tionsofthefreeenergysurfaceforthecorrespondingvari- by a “double-well” structure (see Fig.8), with the wells
ables. Thisisveryusefulbecauseitnotonlyincreasesthe corresponding to the peptide in a “globular” (left mini-
amountofinformationthatcanbegatheredfromagiven mumintheFig.8)andaβ-hairpin(rightminimuminthe
simulation, but it also allows for previously obtained in- Fig.8) folded conformation, respectively. While simple
formationfor a collective variableto be usedto compute enough,the molecule possessesallthe typicalfeatures of
the free energy associated with different variables. For largerpeptide systemsusuallystudiedwithbiomolecular
instance, suppose that in the course of a simulation it simulations. Simulation parameters are as in a previous
becomes apparent that one wishes to address additional study19: the atoms are described by the 1999 version of
questions involving different collective variables. Instead the Cornell et al. force field26, with no cutoff for the
ofstartingfrom“scratch”,onecanre-usethealreadyob- non-bonded interactions. The Berendsen thermostat is
tained biasing potentials and thereby greatly accelerate chosen with τ = 1fs for temperature control. The MD
tp
4
5 5
4 4
4{k
) )
ol ol
m 3{k m
3 3
/ /
al 2{k al 1{k
c c
k k
( ( 2{k
S2 1{k S2
M M 3{k
R R
E E
1 1 4{k
0 0
0 100 200 300 400 0 25 50 75 100 125 150
time (ns) time (ns)
FIG.3: RMSerrorof thefree energyover3.3˚A<Rg <6.3˚A FIG. 4: RMS error of the free energy over 3.3˚A< Rg <6.3˚A
for ABMD simulations at T = 300K with τF = 180 ps (1), formultiplewalkersABMDsimulationswithτF =1nsusingone
360ps(2), 720ps(3) and 1440ps (4). (1), two (2), four (4) and eight (8) trajectories at T =300K.
time step (∆t) is 1 fs for the parallel tempering simula- more accurate than the corresponding MTD run. The
tions, and 2 fs otherwise. ABMDmethodowesitsbetterconvergencetothesmoother
The radius of gyration of the heavy atoms was chosen time evolution of the biasing potential and accurate dis-
to be the collective variable: cretization in Q. The amount of memory used by the
m ABMD simulation to store U values was only ≈ 200×8
Rg = ma (ra−RΣ)2 . (3.1) bytes (considering double pmrecision)for roughly 108 tiny
Xa Σ “hills”thatwereaccumulatedbytheendoftherun. The
Here, R = (m /m )r is the center of mass, with reference MTD simulation with merely 5×103 Gaussians
Σ a a Σ a
mΣ = amPa, and the sum runs over all atoms except required roughly 25 times more memory for the biasing
hydrogePn. The initial configuration is the fully unfolded potential(withonlythepositionsoftheGaussiansstored
peptide. A reference free energy profile, whose error19 explicitly). One canexpect, that ABMDwill be evenmore
in the region of interest is less than 0.15 kcal/mol, was economical when it comes to dealing with multidimen-
computed for benchmarking purposes (see Appendix A sional collective variables, provided that sparse arrays
for details). As a measure of the RMS free energy error, areusedforUm withonlynon-zeroelementsbeingstored
the following construction was used: explicitly.
Although an a priori error estimate for this type of
non-equilibriumsimulationis reallynot feasible, it is ex-
b
ERMS =vuub−1 aZ dξ(cid:16)f1(ξ)−f2(ξ)−∆(cid:17)2, TpehcitsedpotihnattisthileluesrtrroartesdhoinulFdigd.e3c,rweahsiechfosrhoinwcsretahseederτrFor.
ut a for increasing values of τF.
In order to decrease the simulation time required for
where
accuratefreeenergyestimatesevenfurther,themultiple-
b walkervariationof ABMDprovesto be useful. For a mod-
1
∆= dξ f (ξ)−f (ξ) , erate number of walkers, the speedup is nearly linear,
1 2
b−aZ (cid:16) (cid:17) with an additional increase in accuracy coming from the
a
better sampling of the “evolving” canonical distribution
accounts for the arbitrary additive constants in the free (see Fig.4). Parallel tempering improves both the speed
energies f (ξ). Here a = 3.3˚A and b = 6.3˚A, which and the accuracy even more. To this end, we first ran
1,2
correspond to the physical region of interest (see Fig.8). ABMD with τ = 90ps using 2, 4, 6 and 8 replicas at
F
T =300K,331K,365K,403K,445K,492K,543Kand
Figure 2 presents the time dependence of the RMS 600K (during equilibrium MD runs, the peptide configu-
free energy error for the ABMD and reference MTD run19. ration jumps between the two minima on a picosecond
Both simulations have exactly the same kernel width timescale at T = 600K). In all cases, the E was
RMS
4∆ξ = 0.25 ˚A, and flooding timescale τ = 90 ps that found to be ∼ 1 kcal/mol, or less as t → ∞ (data not
F
corresponds to the a posteriori hills acceptance rate re- shown). Again,the improvementinaccuracystems from
ported in Ref.19. It is evident that the AMBD run is the better sampling ofthe “evolving”canonicaldistribu-
5
5 7
-1
-2
4 -3
6
-4
)
ol -5
/m3 1{k ˚A) -6
al ( 5 -7
c g
(k 2{k R -8
S2 -9
M
R 3{k 4 -10
E -11
1 4{k -12
-13
3 -14
0 0 2 4 6 8
0 3 6 9 12 15
NOH kcal/mol
time (ns)
FIG. 6: Free energy map for Ace-GGPGGG-Nme peptide in the
FatIGT.5=: 3R0M0KSefrorrorpoafratlhleelftreeme epneerirnggyoAvBeMrD3s.i3m˚Au<latRiogns<u6s.i3n˚Ag g(saesepRhaefs.e27asfaorfudnecttaiiolsnroefgtahredcinoglletchteivaelvgaorriiathbmlesuRsegdantodNmOakHe.
eight replicas running at T = 300K, 331K, 365K, 403K, this plot).
445K,492K,543K and600K with τF =90ps(1),45ps(2),
22.5ps (3) and 11.25ps(4).
and oxygen atoms and r = 2.5˚A. The sum runs over
0
tion. Then, we ran8 replicaswith smaller values τF and theuniqueO-Hpairs(i.e., eachO-Hpairiscountedonly
were surprised that the accuracy does not degrade, even once), with O and H separated by one or more amino-
for τF =11.2ps as shown in Fig.5. bases along the backbone (27 pairs in total). In other
Finally, we turn to aspects related to the general words,we“re-use”thepreviouslycomputedfreeenergies
replica exchange method and illustrate its potential. As forR togetthefreeenergyinthetwo-dimensionalspace
g
already noted, by using replicas with different collective (R ,N ): the eight first replicas serve as a “sampling
g OH
variables and swapping these at prescribed rates, it is enhancement device” for the ninth replica. The calcula-
possible to obtain projections of the free energy surface tion is carried out in two stages: a “coarse” stage (15ns
for the corresponding variables. It is also possible to with τ = 10ps and 4∆R = 0.25˚A, 4∆N = 0.5)
F g OH
use previously obtained information with respect to one followed by a “fine” stage (50ns with τ = 100ps and
F
collective variable to compute the free energy associated 4∆R = 0.1˚A, 4∆N = 0.25). In both runs exchanges
g OH
with a different variable. For example, suppose that in- between four randomly chosen pairs of replicas were at-
stead of the one-dimensional free energy profile as as a temptedevery100fs. Thefinalfreeenergymapisshown
function of Rg already discussed, one realizes that what in Fig.6. It is clear that this two-dimensional free en-
is actually needed is a two-dimensional profile that in- ergy landscape conveys additional information not con-
cludes information with respect to the number of O-H tained in the one-dimensional free energy plots already
bonds alongthe backbone. The two-dimensionalfree en- discussed. In particular, it allows for a better character-
ergy map is computationally quite expensive, but the ization of the “globular” states of the Ace-GGPGGG-Nme
calculation can be greatly accelerated with the help of : specifically, it is apparent from the Fig.6 that there
the general replica exchange method. We therefore sim- are at least two such states with different values of N
OH
ulated 8+1 = 9 replicas. The eight replicas were run (both correspond to the left minimum in Fig.8). Of
at the previously stated temperatures, with each replica course, one could have re-used the information in the
biased by a static (not evolving) biasing potential corre- one-dimensional R profiles to include other collective
g
sponding to the negated free energy associated with the variables, in addition to N .
OH
radius of gyrationRg at the corresponding temperature, ThegeneralreplicaexchangeABMDmayalsobe advan-
as shown in Fig.8. The additional ninth replica was run
tageous for explicit solvent simulations, which are often
at T = 300K with ABMD flooding in the two collective
notoriously lengthy. Specifically, if one is interested in
variables, i.e., R and the number of O-H bonds along
g the solute and the collective variables do not depend on
the backbone as given by
the solventdegreesoffreedom,then the number ofrepli-
6 cas required to maintain an adequate exchange rate de-
1− r /r
N = OH 0 , pends only very weakly on the amountof solvent(which
OH XO,H1−(cid:0)rOH/r0(cid:1)12 must of course be sufficient as to adequately solvate the
(cid:0) (cid:1) structure), provided that all the replicas are simulated
where r is the distance between a pair of hydrogen at the same temperature. This is because the exchange
OH
6
probability does not explicitly depend on the potential
0
energydifference whenthe temperature ofthe replicasis
the same. While not every choice of collective variables
-3
for different replicas will lead to decent exchange rates, 1{k 2{k
onecanneverthelesstakeadvantageofthispropertyand
-6 3{k
use general replica exchange to enhance the sampling in ol)
a solvated environment. m
/ -9
Inordertodemonstratethe methodinthis regime,we al
c
simulated the Ace-GGPGGG-Nme peptide at T = 300K (k-12
solvated by 171 cyclohexane (C H ) molecules (the to- f
6 12
tal number of atoms was 3139)under periodic boundary -15
conditions using the General28 AMBERForce-Field(GAFF)
for the solvent. We used a truncated octahedron cell of -18
fixedsize(constantvolume)thatcorrespondstotheequi-
3 4 5 6 7
librium density at T = 300K (the equilibrium density R (˚A)
g
value was obtained from a 10ns simulation under con-
stant pressure at T = 300K). The Particle-Mesh Ewald
(PME29)methodwasusedfortheelectrostaticforces,with FIG. 7: Free energy for Ace-GGPGGG-Nme peptide solvated
a 36×36×36 FFT grid and an 8˚A cutoff for the direct in cyclohexane at T = 300K as obtained via: coarse
(non-equilibrium) replica-exchange ABMD (1); finer (non-
sum (same cutoff was used for van der Waals interac- equilibrium) replica-exchange ABMD (2); including the correc-
tions). First, we ran 10 replicas in the “flooding” mode
tion coming from equilibrium biased replica exchange (3).
(i.e., under evolving biasing potentials) using as collec-
tive variables the distances r between the backbone
CC
carbons separated by at least 2 amino-acids (there are
gas phase, the folded β-turn in the cyclohexane solvated
10suchdistancesforAce-GGPGGG-Nme ). Weranfor5ns
peptide is clearly favored over the globular structure.
with τ = 30ps and 4∆r = 1˚A, and then for another
F CC
5nswithτ =150psand4∆r =0.5˚A,attemptingex-
F CC
changesbetweenfiverandomlyselectedpairsevery0.5ps.
IV. CONCLUSIONS AND OUTLOOK.
Asexpected, atthe startofthe simulation,the exchange
rate was nearly 100% decreasing later as the biasing po-
tentials were built up. By the end of the simulation, In summary, we have presented an ABMD method that
when all possible values of the distances had been cov- computesthefreeenergysurfaceofareactioncoordinate
ered, the exchange rate was very disparate between dif- using non-equilibrium dynamics. The method belongs
ferent pairs of replicas. However, for every replica there to the general category of umbrella sampling methods
wasatleastoneotherreplicasuchthattheexchangerate with an evolving potential, and is characterized by only
between them was reasonable(i.e., the whole simulation two control parameters (the flooding timescale and the
did notdegenerateinto ten different non-interactingtra- kernel width) and a favorable O(t) scaling with molec-
jectories). We then set τ =∞ in these 10 replicas, and ular dynamics time t. This scaling can be very impor-
F
added an eleventh replica whose collective variable was tant for large-scale classical MD biomolecular simulations
chosen as the radius of gyration of the heavy atoms. In whenlongsimulationtimesarerequired(see,forexample
otherwords,asbefore,weusethetenreplicasasa“sam- Ref.31, and references therein).
pling enhancementdevice”forthe lastone. We thenran ABMDhasalsobeenextendedtosupportmultiplewalk-
a two-stage flooding scheme: a 5ns coarser stage, with ers and replica exchange. Both variations improve speed
τF = 25ps and 4∆Rg = 0.25˚A; followed by a 10ns finer and accuracy of the method due to the better sampling
stage, with τ = 180ps and 4∆R = 0.2˚A. As before, of the “evolving” canonical distribution. The replica ex-
F g
theexchangeattemptsbetween5randomlyselectedpairs change ABMD has been generalized to include different
ofreplicaswereperformedevery0.5ps. The“raw”ABMD- temperatures and/or collective variables, that move un-
computedfreeenergyassociatedwithR afterthatstage dereitheranevolvingorastaticbiasingpotential. Aside
g
is shown in Fig.7. In a next step, we ran 64 biased sim- from enhancing the sampling, this swapping of replicas
ulations (each comprising of 11 replicas) for 7ns (first hasseveralimportantpracticaladvantages. Mostimpor-
2ns for equilibration followed by 5ns of “production” tantly, it enables one to obtain projections of the free
runs) starting from different initial configurations. We energy surface for any number of collective variables one
set τ =∞ in all replicas (static biasing potentials) and might wish to investigate. In addition, one can re-use
F
recorded the values of R in the eleventh replica every previously obtained results in order to enhance the sam-
g
10picoseconds. Wethenusedthelog-splinealgorithmof pling of new collective variables. It is also possible to
Stone30 at. al. toestimatethe(biased)log-densityofthe exploit the fact that exchangerates at the same temper-
R valuesatequilibrium. Thisledustothefinalshapeof atureareindependentofthepotentialenergytoenhance
g
the free energy curve shown in Fig.7. Compared to the samplingofasoluteinaminimumamountofsolvent(for
7
addition we thank NC State HPC for computational re-
-16 sources.
-18
)
ol-20
m APPENDIX A: REFERENCE FREE ENERGY
/ CURVE.
al
c-22
k
(
f Here,weprovidesimulationdetailswithregardstothe
-24
reference free energy curve. We first ran short (5 ns,
eightwalkerswith τ =60psand4∆ξ =0.2˚A)multiple-
F
-26 walkersABMDatT =600Ktoreconstructtheglobalwell.
ThiswasfollowedbyparalleltemperingABMDruns,using
3.25 3.5 3.75 4 4.25 4.5 4.75 5 5.25
R (˚A) the biasing potential obtained from the multiple-walkers
g
simulation as the zero-time value for the biasing poten-
tials at different temperatures. We used eight replicas
FIG. 8: The accurate free energies of Ace-GGPGGG-Nme pep-
at T = 300K, 331K, 365K, 403K, 445K, 492K, 543K
tideingasphaseasfunctionofRgatT =300K,331K,365K, and 600K and attempted exchanges every 100 MD steps
403K, 445K,492K, 543K and 600K (from bottom to top).
(0.1ps). The simulation started with τ = 60ps and
F
4∆ξ = 0.2˚A, and ran for 105 exchanges. We then set
collective variables independent of solvent atom coordi- τF to 600ps, 4∆ξ to 0.1˚A and ran for 5×105 more ex-
nates). We have implemented the ABMD method in the changes. Finally, this was followedup with 1×106 more
AMBERpackage9,andplantodistributeitfreely. Here,we exchanges with τF =6ns and 4∆ξ =0.1˚A.
havedemonstratedtheworkingsoftheABMDmethodwith We then ran a very long (3 × 107 exchanges, 0.1ps
a study of the folding of the Ace-GGPGGG-Nme peptide, between exchanges) biased parallel tempering simula-
The application of ABMD to more complicated biomolec- tion in the spirit of Ref.19, in order to get an a pos-
ular systems is reserved for future publications. teriori error estimate. From the resulting histogram it
follows that the error does not exceed ≈ 0.15kcal/mol
for 3.3˚A< R <6.3˚A. The RMS error is probably much
g
Acknowledgments smaller, since 0.15 corresponds to the absolute non-
uniformity of the histogram, i.e., the maximum error,
This research was partly supported by NSF under over 3.3˚A< R <6.3˚A. The accurate free energy curves
g
grants ITR-0121361 and CAREER DMR-0348039. In as a function of temperature are shown in Fig.8.
∗ Electronic address: [email protected] 10 B.Ensing,A.Laio,F.L.Gervasio,M.Parrinello,andM.L.
1 D. Frenkel and B. Smit, Understanding Molecular Sim- Klein, J. Am.Chem. Soc. 126, 9492 (2004).
ulation, Computational Science Series (Academic Press, 11 S. V. Churakov, M. Ianuzzi, and M. Parrinello, J. Phys.
2002). Chem. B 108, 11567 (2004).
2 T.Huber,A.E.Torda,andW.F.vanGunsteren,J.Com- 12 F. Gervasio, A. Laio, and M. Parrinello, J. Am. Chem.
put.Aided. Mol. Des. 8, 695 (1994). Soc. 124, 2600 (2005).
3 F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 13 M.Ceccarelli,C.Danelon,A.Laio,andM.Parrinello,Bio-
(2001). phys.J. 87, 58 (2004).
4 E. Darve and A. Pohorille, J. Chem. Phys. 115, 9169 14 M.IannuzziandM.Parrinello,Phys.Rev.Lett.93,025901
(2001). (2004).
5 A.LaioandM.Parrinello,Proc.Natl.Acad.Sci.99,12562 15 A. L. A. Stirling, M. Ianuzzi and M. Parrinello,
(2002). ChemPhysChem 5, 1558 (2004).
6 M. Iannuzzi, A. Laio, and M. Parrinello, Phys. Rev. Lett. 16 E. Asciutto and C. Sagui, J. Phys. Chem. A 109, 7682
90, 238302 (2003). (2005).
7 T. Leli`evre, M. Rousset, and G. Stoltz, J. Chem. Phys. 17 J. G. Lee, E. Asciutto, V. Babin, C. Sagui, T. A.Darden,
126, 134111 (2007). and C. Roland, J. Phys. Chem. B 110, 2325 (2006).
8 G.Bussi,A.Laio, andM.Parrinello, Phys.Rev.Lett. 96, 18 T.Ikeda,M.Hirata,andT.Kimura,J.Chem.Phys.122,
090601 (2006). 244507 (2005).
9 D. A. Case, T. E. Cheatham, III, T. Darden, H. Gohlke, 19 V.Babin,C.Roland,T.A.Darden,andC.Sagui,J.Chem.
R. Luo, K. M. Merz, Jr., A. Onufriev, C. Simmerling, Phys. 125, 204909 (2006).
B. Wang, and R. Woods, J. Computat. Chem. 26, 1668 20 A. Laio, A. Rodriguez-Fortea, F. L. Gervasio, M. Cec-
(2005). carelli, and M. Parrinello, J. Phys. Chem. B 109, 6714
8
(2005). Caldwell, and P. A. Kollman, J. Am. Chem. Soc. 117,
21 B. W. Silverman, Density Estimation for Statistics and 5179 (1995).
DataAnalysis,Monographsonstatisticsandappliedprob- 27 A.Preusser,ACMTransactionsonMathematicalSoftware
ability (Chapman and Hall, 1986). 15, 79 (1989).
22 P. Raiteri, A. Laio, F. L. Gervasio, C. Micheletti, and 28 J. Wang, R. Wolf, J. Caldwell, P. Kollman, and D. Case,
M. Parrinello, J. Phys.Chem. 110, 3533 (2006). J. Comp. Chem. 25, 1157 (2004).
23 Y.Sugita,A.Kitao,andY.Okamoto,J.Chem.Phys.113, 29 T. A. Darden, D. M. York, and L. G. Pedersen, J. Chem.
6042 (2000). Phys. 98, 10089 (1993).
24 G. Bussi, F. L. Gervasio, A. Laio, and M. Parrinello, J. 30 C.J.Stone,M.Hansen,C.Kooperberg,andY.K.Truong,
Am. Chem. Soc. 128, 13435 (2006). Annals of Statistics 25, 1371 (1997).
25 S.PianaandA.Laio,J.Phys.Chem.B111,4553(2007). 31 M.R.Shirts,J.W.Pitera,W.C.Swope,andV.S.Pande,
26 W.D.Cornell, P.Cieplak, C. I.Bayly,I.R.Gould, K.M. J. Chem. Phys.119, 5740 (2003).
Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W.