Annealed importance sampling of dileucine peptide
7 Annealed importance sampling of dileucine
peptide
Edward Lyman∗, and Daniel M. Zuckerman†
Department of Computational Biology, School of Medicine
3079 BST3 , 3501 Fifth Ave.,
University of Pittsburgh, Pittsburgh, PA 15261
October 28, 2018
Abstract
Annealed importance sampling is a means to assign equilibrium
weights to a nonequilibrium sample that was generated by a simu-
lated annealing protocol[1]. The weights may then be used to calculate
equilibrium averages, and also serve as an “adiabatic signature” of the
chosen cooling schedule. In this paper we demonstrate the method on
the 50-atom dileucine peptide, showing that equilibrium distributions
are attained for manageable cooling schedules. For this system, as
∗elyman@ccbb.pitt.edu
†dmz@ccbb.pitt.edu
http://arxiv.org/abs/0704.0464v1
näıvely implemented here, the method is modestly more efficient than
constant temperature simulation. However, the method is worth con-
sidering whenever any simulated heating or cooling is performed (as is
often done at the beginning of a simulation project, or during an NMR
structure calculation), as it is simple to implement and requires mini-
mal additional CPU expense. Furthermore, the näıve implementation
presented here can be improved.
1 Introduction
Simulated annealing (SA) is used in a wide variety of biomolecular
calculations. Crystallographic refinement protocols[2] and standard
NMR structure calculations[3, 4, 5] both rely on SA to optimize a
“target function,” constructed so that the global minimum of the tar-
get function corresponds to the native structure. Molecular dynamics
calculations often begin by cooling a configuration from a high tem-
perature ensemble to a lower temperature, at which the simulation is
to be performed.
In this paper, we consider a different use for SA calculations. Since
a set of structures that is generated by a series of SA trajectories is a
nonequilibrium sample, they may not be used to calculate equilibrium
averages. However, Neal demonstrated a simple procedure, called “an-
nealed importance sampling” (AIS) that allows the nonequilibrium
sample to be reweighted into an equilibrium one[1]. AIS is closely
connected with the Jarzynski relation[6]. To our knowledge, the al-
gorithm has only appeared once in the chemical physics literature[7],
where it was used (along with sophisticated Monte Carlo techniques)
to sample a one-dimensional potential. Here, we demonstrate an ap-
plication of the AIS algorithm to generate an equilibrium sample of
an implicitly solvated peptide, and discuss other uses for AIS which
may of interest to the molecular simulation community.
The basic idea which underlies SA is also the motivation for other
temperature based sampling methods, notably J-walking[8], simulated
tempering[9, 10] and replica exchange/parallel tempering[11, 12]. By
coupling a simulation to a high temperature reservoir, it is hoped that
the low temperature simulation may explore the configuration space
more thoroughly. This is achieved by thermally activated crossing
of energetic barriers, which are large compared to the thermal en-
ergy scale of the lower temperature simulation, but are crossed more
frequently at higher temperature. Simulated and parallel tempering
differ in the way that the different temperature simulations are cou-
pled. Simulated tempering heats and then cools the system, in a way
that maintains an equilibrium distribution. Parallel tempering couples
simulations run in parallel at different temperatures by occasionally
swapping configuartions between temperatures, again in such a way
that canonical sampling is maintained.
AIS offers yet another approach to utilizing a high temperature
ensemble for equilibrium sampling at a lower temperature. A sam-
ple of a high temperature ensemble is annealed to a lower temper-
ature, by alternating constant temperature simulation with steps in
which the tempertaure is jumped to a lower value. Each annealed
structure is assigned a weight, which depends on the trajectory that
was traced during the annealing process. Equilibrium averages over
the lower temperature ensemble may then be calculated by a simple
weighted average. Furthermore, the distribution of trajectory weights
contains useful information about the statistics of the annealed sam-
ple. Roughly, a schedule which quenches high temperature structures
very rapidly to low temperature will result in a sample dominated by
a few high weight structures, resulting in poor statistics. This con-
nection between the distribution of weights and the extent to which
the schedule is not adiabatic ought to be of interest to anyone who
uses SA protocols—whether for equilibrium sampling or for structure
calculation.
We have used the AIS method to generate 298 K equilibrium en-
sembles of the dileucine peptide, by annealing structures from a 500
K distribution with several different cooling schedules. For the most
efficient schedule used, we found a modest gain (about a factor of 3)
over constant temperature simulation. This result is consistent with
earlier observations on the expected efficiency of temperature-based
sampling methods[13].
2 Theory
Consider a standard simulated annealing (SA) trajectory, in which a
protein is slowly cooled from a conformation x at a (high) temperature
TM . The cooling is achieved by alternating constant temperature
dynamics with “temperature jumps,” during which the temperature
is lowered instantaneously. Usually, the system is cooled to a low
temperature, since the aim of standard SA calculations is to find the
global minimum on the energy landscape. But we can imagine instead
ending the run at T0 = 300 K—in fact, we can think of many such
runs, all ending at 300 K. We then have an ensemble of conformations,
though clearly not distributed canonically at T0. We would like to
know if there is a way to reweight this distribution, so that it can be
used to compute equilibrium averages at T0. The affirmative answer
is provided by the annealed importance sampling (AIS) method.
To make the discussion more concrete, consider many independent
annealing trajectories xj(t) which at time tM−1 have just been cooled
from inverse temperature βM to βM−1. As usual, each temperature
defines a distribution of conformations: πi(x) ∝ exp[−βiU(x)]. Imme-
diately after tM−1, before the system is allowed to relax to πM−1(x),
we can compute the equilibrium average of an arbitrary quantity A
over πM−1(x) by using the weight w(x) = πM−1(x)/πM (x):
〈A〉M−1ZM−1 =
dxA(x)πM−1(x) =
dxA(x)πM (x)w(x), (1)
where 〈A〉i denotes an average over πi, and Zi =
dxπi(x). In other
words, we may reweight the distribution πM(x) to calculate averages
over πM−1(x), by multiplying by the ratio of Boltzmann factors.
Generalizing the argument toM temperature steps is straightforward[1],
by forming the product of weights for successive cooling steps:
wj ≡ w(xj(t0)) =
πi−1(xj(ti−1))
πi(xj(ti−1))
. (2)
Equation 2 gives the weight for trajectory j, cooled at successive times
tM−1, tM−2,... through inverse temperatures βM , βM−1,... to reach
conformation xj(t0). At each temperature, reweighting ensures that
averages may be calculated for the appropriate canonical distribution,
even though the system has not yet relaxed.
The AIS idea is easily turned into an algorithm for producing a
canonical distribution from serially generated annealing trajectories:
(i) Generate a sample of the distribution πM(x), by a
sufficiently long simulation at TM .
(ii) Pull a conformation from πM (x) at random and anneal
down to β0, yielding conformation x1(t0). Keep track of the
weight w(x1(t0)) for this trajectory by Eq. 2.
(iii) Repeat steps (iii) and (iv) N times, yielding congiura-
tions xj and weights w(xj) ≡ wj for j = 1, 1, ..., N .
Equilibrium averages at temperature T0 are then calculated by a
weighted average:
〈A〉0 =
j=1wjAj∑N
The cooling schedule is defined by the number and spacing of the
temperature steps, as well as the duration of the constant temperature
simulation at each step. As available resources necessarily limit the
CPU time spent on each annealing trajectory, careful consideration of
the schedule is in order. Clearly, a schedule in which high temperature
configurations are quenched in one step to low temperature amounts
to a single-step reweighting procedure[14]. We may expect that such a
schedule would be quite ineffective for large temperature jumps, since
very few configurations in the high temperature distribution have ap-
preciable weight in the low temperature distribution. By introducing
intermediate steps, the system is allowed to relax locally, bridging the
high and low temperature distributions in a way that echoes replica
exchange protocols[11, 12], simulated tempering[9, 10], and the multi-
ple histogram method[15]. However, the “top-down” structure of the
algorithm most closely resembles J-walking[8, 16].
3 Results
The dileucine peptide (ACE-[Leu]2-NME) is good choice for the vali-
dation of new algorithms, as it is small enough (50 atoms, including
nonpolar hydrogens) that exhaustive sampling by standard simulation
methods is possible, yet more akin to protein systems than a one- or
two-dimensional “toy” model.
The high temperature ensemble was generated by 300 nsec of
Langevin dynamics at TM = 500 K, as implemented in Tinker v.
4.2.2[17], with a timestep of 1.0 fsec, and a friction constant of 91
psec−1, and solvation was treated by the GB/SA method[18]. Frames
were written every psec, resulting in a sample of 3× 104 frames in the
high temperature sample.
The 500 K sample was annealed down to 298 K using 4 different
schedules, consisting of a total of 3, 5, 9, and 17 temperatures, includ-
ing the endpoints. In each case, the temperatures were distributed
geometrically. Following each temperature jump, the velocities were
reinitialized by sampling randomly from the Maxwell-Boltzmann dis-
tribution, and then allowed to relax at constant temperature for a
time tR = 0.5 psec (except where noted) with the protocol described
above. A total of N = 1.6× 104 annealing trajectories were generated
for each schedule. The control of the integration routine to effect the
annealing, as well as the calculation of the trajectory weights, were
implemented in a Perl script.
Figure 1 shows that the 298 K distribution of energy is recovered
by the AIS procedure. It is noteworthy that the 500 K distribution
(corresponding to the high T sample) overlaps very little with the 298
K distribution, and yet the 298 K distribution is reproduced well for
the two slowest schedules. Equally interesting is how poorly the algo-
rithm performs when the structures are cooled too rapidly, especially
on the low E side of the distribution, where there is no overlap with
the high T distribution. We conclude that the schedules with 3 or 5
T -steps quench the structures too rapidly, resulting in many of the tra-
jectories becoming “stuck” in high-energy states that are metastable
at 298 K.
This last observation may be quantified by asking, “How many
of the annealed structures contribute appreciable weight to averages
calculated with Eq. 3?” To address this question, for each schedule we
estimated the number of configurations n which contribute appreciable
weight to the averages:
≡ fN, (4)
where wmax is the largest weight observed (see Table 1). If this number
is near 1, then a small number of trajectories dominate the average—
see Eq. 3 —and poor results should be expected. The effective fraction
of the annealing trajectories which generate “useful” or “successful”
structures is denoteed by f .
A more complete picture is provided by the full distribution of
the (logarithm of) trajectory weights (Fig. 2). For each schedule, the
weights which contribute the most to the T = 298 K sample are to
the right, at large values of w. The trend is clear—as slower cooling
is effected, the distribution narrows and shifts to the right. It has
been shown that the accuracy of averages computed from this type of
protocol is roughly related to the variance of the (adjusted) weights[1].
(The adjusted weight is the weight divided by the average weight.)
This “rule of thumb” is borne out by the data in Fig. 2 and Table
1—as the cooling slows down the distribution of weights narrows, and
the number of trajectories contributing to the equilibrium averages
increases. This type of analysis may serve as a means of distinguishing
between annealing schedules to decide on a cooling schedule which is
slow enough to yield reasonable estimates of equilibrium averages. It
is also essential for optimizing an AIS protocol for sampling efficiency,
as discussed in the next few paragraphs.
How much better than standard simulation (if at all) is equilib-
rium sampling by AIS? In order to make a direct comparison between
AIS and constant temperature simulation, we need to compare the
CPU time invested per statistically independent configuration in each
protocol. For the constant temperature simulation, this time may be
estimated in several ways[19, 20], and is essentially the time needed
for the simulation to “forget” where it has been. Following the con-
vention for correlation times, we call this time τi = τ(Ti), where i
labels the temperature: M for the high T distribution, and 0 for the
low T distribution. For the system studied here, τM = 0.8 nsec and
τ0 = 3.0 nsec, as estimated from timseries of the α → β backbone
dihedral transition[13].
The total cost to generate a structure in an AIS simulation is the
sum of the costs of generating a structure in the high T distribution
plus that for the annealing phase. Of course, not every annealing
trajectory contributes to thermodynamic averages(Eq. 3). What then
is the total cost tcost of a “successful” annealed structure? The first
part is from high temperature sampling—i.e., τM . The second part is
the cost of all the annealing trajectories, divided by the number which
contribute to equilibrium averages. The time tanneal is the time spent
annealing each structure:
tanneal = tR(M − 2) (5)
Recall that tR is the duration of the constant temperature relaxation
steps, and there is no relaxation phase at the highest and lowest tem-
peratures.
The total cost tcost is then the sum of τM and tanneal:
tc = τM + tanneal/f. (6)
The efficiency of an AIS protocol may then be computed by taking
the ratio R ≡ τ0/tcost (see Table 1), which gives the factor by which
an AIS protocol is more or less efficient than constant temperature
simulation. The data in Table 1 show that the best schedule used
here offer a modest speedup over constant temperature simulation, of a
factor of about 3. These findings are in agreement with an analysis we
have published of another temperature-based sampling protocol[13].
We note that an optimized AIS protocol would require tuningN based
on (perhaps preliminary) estimates of f .
It is instructive to compare the AIS results to simple reweighting—
i.e., AIS with no intermediate temperature steps or relaxation. In
this case, no computer time is spent annealing, and the efficiency
gain is simply τ0/τM = 3.75. The fraction f is of course reduced
compared to any AIS protocol—when reweighting our 500 K dileucine
trajectory to 298 K distribution, f = 1.3 × 10−4—but this has no
impact on the efficiency, provided a sufficient number of snapshots are
available for reweighting. However, it is clear that f will be greatly
reduced in systems which undergo a folding transition upon lowering
the temperature. This is simply a reflection of the fact that there is
negligible overlap between the folded and unfolded distributions. In
such cases, a useful reweighting protocol would require the generation
of astronomical numbers of structures in the TM distribution, and
annealing is advised.
4 Conclusion
We have demonstrated the application of Neal’s annealed importance
sampling (AIS) algorithm for equilibrium sampling of the dileucine
peptide. AIS allows the calculation of equilibrium averages from a
nonequilibrium sample of strutures that results from a simulated an-
nealing protocol. To our knowledge, AIS has not previously been
applied to a molecular system. While the method, as näıvely im-
plemented here, represents only a modest improvement over constant
temperature simulation, it is interesting for several reasons beyond
equilibrium sampling.
First, in applications where simulated annealing is already in widespread
use (most notably, NMR structure calculations[3, 4, 5]), the path
weights may be used to calculate (perhaps noisy) equilibrium aver-
ages, and perhaps ultimately Boltzmann-distributed ensembles. The
path weights also contain information that can be used to discriminate
between different schedules, which may provide a way to optimize the
schedule, based on the analysis of tann, the cost of annealing to “good”
structures.
Second, it may be possible to improve considerably on the efficiency
of the method by implementing a more sophisticated version, which
uses a resampling procedure to prune the low weight paths at each
cooling step. (For a detailed discussion of resampling methods, see
the book by Liu[21].) In this approach, we first cool some number N
of structures from the high temperature (TM ) ensemble, yielding N
weighted structures at TM−1. We then resample N times from this
TM−1 ensemble, according to the cumulative distribution function of
the weights, pruning the low weight paths without biasing the sample.
This type of approach was recently applied successfully to sampling
near native protein configurations of a discretized and coarse-grained
model[22]. Nevertheless, we emphasize that the ultimate efficiency of
any AIS protocol limited by the intrinsic sampling rate of the highest
temperature, which may be modest; see Ref. ??.
Finally, the AIS procedure could be naturally combined with “an-
nealing” in the parameters of the Hamiltonian. Such a hybrid of AIS
and Hamiltonian switching might be used, for example, to transform
an NMR target function into a molecular mechanics potential function,
over the course of a structure calculation. The result of such a cal-
culation would be an equilibrium ensemble of structures, distributed
according to the molecular mechanics potential. Such ensembles would
find wide application, for instance in docking or homology modeling.
Acknowledgements The authors thank Gordon Rule for several
enlightening discussions about NMR methodology. D. Z. thanks Chris
Jarzynski for alerting him to Neal’s work on AIS. This research was
supported by the NSF (MCB-0643456), the NIH (GM076569), and
the Department of Computational Biology, University of Pittsburgh.
References
[1] Radford M. Neal. Annealed importance sampling. Stat. and
Comp., 11:125–139, 2001.
[2] Axel T. Brünger, Paul D. Adams, G. Marius Clore, Warren L.
DeLano, Piet Gros, Ralf W. Grosse-Kuntsleve, Jian-Sheng Jiang,
John Kuszewski, Michael Nilges, Navraj S. Pannu, Randy J.
Read, Luke M. Rice, Thomas Simonson, and Gregory L. War-
ren. Crystallograhy and NMR system: a new software suite for
macromolecular structure determination. Acta. Cryst., D54:905–
921, 1998.
[3] C.D. Schwieters, J.J. Kuszewski, N. Tjandra, and G.M. Clore.
The Xplor-NIH NMR molecular structure determination package.
J. Magn. Res., 160:66–74, 2003.
[4] P. Güntert, W. Braun, and K. Wüthrich. Torsion angle dynamics
for NMR structure calculation with the new program DYANA.
J. Mol. Biol., 273:283–298, 1997.
[5] Axel T. Brünger, Paul D. Adams, and Luke M. Rice. New ap-
plications of simulated annealing in X-ray crystallography and
solution NMR. Structure, 15:325–336, 1997.
[6] C. Jarzynski. Nonequilibrium equality for free energy differences.
Phys. Rev. Lett., 78:2690–2693, 1997.
[7] S. Brown and T. Head-Gordon. Cool walking: A new Markov
chain Monte Carlo method. J. Comp. Chem., 24:68–76, 2002.
[8] D. D. Frantz, D. L. Freeman, and J. D. Doll. Reducing quasi-
ergodic behavior in Monte Carlo simulations by J-walking: appli-
cations to atomic clusters. J. Chem. Phys., 93:2768–2783, 1990.
[9] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N.
Vorontsov-Velyaminov. New approach to Monte Carlo calculation
of the free energy: Method of expanded ensembles. J. Chem.
Phys., 96:1776–1783, 1992.
[10] E. Marinari and G. Parisi. Europhys. Lett., 19:451–458, 1992.
[11] Charles J. Geyer. Markov chain Monte Carlo maximum likeli-
hood. In E. M. Keramidas, editor, Proceedings of the 23rd sympo-
sium on the interface, Computing science and statistics. Interface
foundation of North America, 1991.
[12] David J. Earl and Michael W. Deem. Parallel tempering: theory,
applications, and new perspectives. Phys. Chem. Chem. Phys.,
23:3910–3916, 2005.
[13] Daniel M. Zuckerman and Edward Lyman. A second look at
canonical sampling of biomolecules using replica exchange simu-
lation. J. Chem. Th. and Comp., 4:1200–1202, 2006.
[14] Alan M. Ferrenberg and Robert H. Swendsen. New Monte
Carlo technique for studying phase transitions. Phys. Rev. Lett.,
61:2635–2638, 1988.
[15] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and
J. M. Rosenberg. The weighted histogram analysis method for
free energy calculations in biomolecules. i. the method. J. Com-
put. Chem., 13:1011–1021, 1992.
[16] Alexander Matro, David L. Freeman, and Robert Q. Topper.
Computational study of the structures and thermodynamic prop-
erties of ammonium chloride clusters using a parallel jump-
walking approach. J. Chem. Phys., 104, 1996.
[17] http://dasher.wustl.edu/tinker/.
[18] W. C. Still, A. Tempczyk, and R. C. Hawley. Semianalytical
treatment of solvation for molecular mechanics and dynamics. J.
Am. Chem. Soc., 112:6127–6129, 1990.
[19] A. M. Ferrenberg, D. P. Landau, and K. Binder. Statistical and
systematic errors in monte carlo sampling. J. Stat. Phys., 63:867–
882, 1991.
[20] Edward Lyman and Daniel M. Zuckerman. On the convergence
of biomolecular simulations by evaluation of the effective sample
size. preprint: http://xxx.lanl.gov/abs/q-bio.QM/0607037.
[21] Jun S. Liu. Monte Carlo strategies in scientific computing.
Springer, New York, 2001.
[22] Jinfeng Zhang, Ming Lin, Rong Chen, Jie Liang, and Jun S. Liu.
Monte Carlo sampling of near-native structures of proteins with
applications. PROTEINS, 66:61–68, 2007.
T-steps Annealing time Successful Fractional Net cost Efficiency
tanneal structures success rate gain
M = (M − 2)tR n f ≡ n/N tcost (nsec) R
3† 0.5 psec 7.1 4.4× 10−4 1.94 1.5
5† 1.5 psec 43.7 2.7× 10−3 1.36 2.2
17† 1.5 psec 137.6 8.6× 10−3 0.97 3.1
33 1.5 psec 46.2 2.9× 10−3 1.32 2.3
9† 3.5 psec 163.2 1.0× 10−2 1.15 2.6
17 7.5 psec 205.3 1.3× 10−2 1.38 2.2
17 15.0 psec 237.2 1.5× 10−2 1.80 1.7
17 30.0 psec 353.8 2.2× 10−2 2.16 1.4
Table 1: Comparison of the efficiency of AIS between several cooling sched-
ules. n is given by Eq. 4, tcost is given by Eq. 6. The efficiency gain is the total
simulation time invested in each successful annealed structure tcost divided
by the time needed to generate an indepenendent structure by constant tem-
perature simulation. The † indicates schedules for which data are presented
in Figs. 1 and 2.
Figure Legends
Figure 1.
Distribution of energies, from standard, constant temperature simu-
lation and AIS. The dashed line is the T = 500 K distribution that
was used for the high T ensemble. The other data compare a 300
nsec, T = 298 K constant temperature simulation to 298 K ensembles
generated by the AIS algorithm with different cooling schedules. The
schedules are discussed in Table 1.
Figure 2.
Distribution of the logarithm of trajectory weights for the four cooling
schedules used in Fig. 1 and discussed in Table 1.
Figure 1:
Figure 2:
	Introduction
	Theory
	Results
	Conclusion
