Distributed Replica Sampling
Tomas Rodinger,†,‡,§ P. Lynne Howell,†,‡ and Re´gis Pome`s*,†,‡
Structural Biology and Biochemistry, The Hospital for Sick Children, 555 UniVersity
AVe., Toronto, Ontario, Canada M5G 1X8, Department of Biochemistry, UniVersity of
Toronto, 1 King’s College Circle, Toronto, Ontario, Canada M5S 1A8, and Institute of
Biomaterials and Biomedical Engineering, UniVersity of Toronto, 4 Taddle Creek
Road, Rm. 407, Rosebrugh Building, Toronto, Ontario, Canada M5S 3G9
Received December 1, 2005
Abstract: We present a simple and general scheme for efficient Boltzmann sampling of
conformational space by computer simulation. Multiple replicas of the system differing in
temperature T or reaction coordinate λ are simulated independently. In addition, occasional
stochastic moves of individual replicas in T or λ space are considered one at a time on the
basis of a generalized Hamiltonian containing an extra potential energy term or bias that depends
on the distribution of all replicas. The algorithm is inherently suited for shared or heterogeneous
computing platforms such as a distributed network.
Introduction
Despite rapid advances in computer technology, statistical
sampling of systems governed by rugged potential energy
surfaces remains a challenge, which underscores the need
for efficient computational algorithms.1-6 For example, the
time scales of conformational transitions of biomolecular
systems span many orders of magnitude. Many of these
events are not directly accessible by atomistic computer
simulations, precluding Boltzmann sampling of conforma-
tional space.
A number of statistical mechanical methods using general-
ized ensembles have been developed in recent years to
enhance sampling efficiency. Some of these approaches are
designed to induce a random walk either in temperature,
potential energy, spatial coordinates, or transformation
parameters in the Hamiltonian. In simulated tempering
(ST),7,8 a random walk in temperature is achieved by taking
periodic steps from absolute inverse temperature "i to "j
subject to acceptance of the following Metropolis Monte
Carlo test:
where "i and "j are the inverse of the Boltzmann constant
multiplied by absolute temperature, the potential energy E(q)
is a function of atomic coordinates q, and ai and aj are weight
factors that must be determined in advance using an iterative
procedure. Similarly, the multicanonical algorithm (MU-
CA)9,10 results in a random walk in potential energy space.
Extending the Hamiltonian to include a transformation
coordinate λ is another way to get around barriers because
the higher dimensionality of phase space results in an
increased number of alternate routes by which barriers can
be avoided. A random walk along a transformation coordi-
nate can be realized via adaptive umbrella sampling (AUS)
methods,11-15 where the Hamiltonian governing movement
is defined as
In eq 2, the potential energy E is a function of atomic
coordinates q and transformation coordinate λ, and U(λ) is
a biasing potential that must be adapted such that it removes
any barriers that exist along the transformation coordinate.
In ST, MUCA, or AUS, random walk behavior is achieved
only once the system “adapts” or learns the correct weight
factors for movement along the parameter of interest, a
system-dependent procedure that introduces extra complexity
as well as possible artifacts if not done carefully (see ref 16
* Corresponding author tel.: 416-813-5686; fax: 416-813-5022;
e-mail: pomes@sickkids.ca.
† The Hospital for Sick Children.
‡ Department of Biochemistry, University of Toronto.
§ Institute of Biomaterials and Biomedical Engineering, University
of Toronto.
paccept ) min(1, exp{-[("j - "i)E(q) + (aj - ai)]}) (1)
H(q,λ) ) E(q,λ) + U(λ) (2)
725J. Chem. Theory Comput. 2006, 2, 725-731
10.1021/ct050302x CCC: $33.50 © 2006 American Chemical Society
Published on Web 03/22/2006

for more details). Furthermore, although adaptive methods
clearly improve sampling uniformity along the coordinate,
in some cases, uniform sampling may only be achieved in
very long simulations because of ruggedness in orthogonal
degrees of freedom.17
In the past decade, replica exchange (RE),18-20 also known
as multiple Markov chain or parallel tempering, has emerged
as a powerful alternative to adaptive methods. In RE, multiple
noninteracting copies or replicas of a system, each governed
by the same potential energy function but differing in
temperature, are simulated at once. Periodically, the simula-
tions are halted, and replicas i and j with neighboring
temperatures Ti and Tj are swapped with a probability given
by the following Metropolis Monte Carlo condition:
Temperature swaps cause each replica to undergo a random
walk in temperature, providing a means (for any replica) to
escape from a potential energy trap. Like adaptive methods,
RE may also be used to attain a random walk in a parameter
of the Hamiltonian,21-23 λ. Swaps between adjacent replicas
i and j are accepted with a probability given by
When data from all replicas are taken together, perfect
sampling uniformity along the temperature or reaction
coordinate is attained. Furthermore, no adaptation of a biasing
potential is required because the weight factors are known a
priori. In the realm of biomolecular studies, large-scale
applications of RE have pushed back limitations in the scope
of problems and the size of systems studied.19,23-31 Recent
studies have applied the method to folding a 46 amino acid
protein domain,32 and the approach has been extended to
two-dimensional random walks in pressure and temperature.31
The principal drawback of RE is that it requires a homoge-
neous and often large computing cluster to implement
efficiently. Because the number of replicas required follows
the square root of the number of degrees of freedom in the
system,18 the simulation of complex systems eventually
becomes impractical.
Here, we present distributed replica (DR) sampling, a
simple and general scheme for efficient Boltzmann sampling
of conformational space by computer simulations. As in RE,
multiple replicas of the system, covering a preassigned range
in temperature or reaction coordinate space, are simulated
independently. The main difference with RE is that, instead
of performing pairwise exchanges of replicas, stochastic
moves of individual replicas in T or λ space are considered
one at a time. The target distribution of replicas is enforced
by a generalized Hamiltonian containing an extra potential
energy bias that depends on the distribution of all replicas.
Contrary to adaptive methods, DR does not require prelimi-
nary determination of the bias. We show with a simple but
relevant example that the DR algorithm leads to random-
walk movement along the parameter of interest with an
efficiency comparable to that of RE. Because it avoids the
need for all replicas to run synchronously, DR is intrinsically
suited for implementation on a shared or inhomogeneous
computing platform such as a large-scale distributed network.
Method
Let us consider N noninteracting copies (or “replicas”) of a
system of interest, all of which are governed by an identical
potential energy function, E(qi,λi), where qi represents the
atomic coordinates of the atoms in replica i and λi is the
coupling parameter for the reaction coordinate of interest
(the reaction in question may be either an alchemical or a
spatial transformation). The DR method makes use of an
additional potential energy term, D(λ1,λ2,...,λΝ), henceforth
referred to as the distributed replica potential energy (DRPE),
which enforces the distribution of replicas across the range
of the transformation coordinate (i.e., an energy penalty is
associated with a nonideal distribution). The generalized
Hamiltonian for state X ) {q1,λ1,q2,λ2,...,qN,λN}, combining
all replicas together with the DRPE, is given by
The weight factor for state X is given by
We consider one λ move at a time. Suppose that the λ value
of replica m is to be changed from λm to λm + δλm, thus
taking state X to state X :
In order for the exchange process to converge toward the
equilibrium distribution, it is sufficient to impose the detailed
balance condition on the transition probability p(XfX ):
From eqs 5, 6, and 8, we have
where
This can be satisfied using the Metropolis Monte Carlo
criterion:
Like RE, the DR method may also be performed in
temperature space. To this end, we first define X )
{q1,"1,...,qN,"N}. The corresponding dimensionless general-
ized Hamiltonian is
Hλ(X) )
i)1
N
E(qi,λi) + D(λ1,λ2,...,λN) (5)
W(X) ) exp[-" Hλ(X)] (6)
X ) {q1,λ1,...,qm,λm,...,qN,λN}f X )
{q1,λ1,...,qm,λm + δλm,...,qN,λN} (7)
W(X) p(XfX ) ) W(X ) p(XfX) (8)
p(XfX )
p(XfX) ) exp(-"∆) (9)
∆ ) E(qm,λm + δλm) - E(qm,λm) + D(λ1,λ2,...,λm +
δλm,...,λN) - D(λ1,λ2,...,λm,...,λN) (10)
paccept ) min[1, exp(-"∆)] (11)
H"
*(X) )
i)1
N
"i E(qi) + D*("1,"2,...,"N) (12)
paccept ) min(1, exp{-["i - "j][E(qj) - E(qi)]}) (3)
paccept ) min(1, exp{-"[E(qi,λj) - E(qi,λi) + E(qj,λi) -
E(qj,λj)]}) (4)
726 J. Chem. Theory Comput., Vol. 2, No. 3, 2006 Rodinger et al.