Zhongwei Zhu and Mark E. Tuckerman
Dept. of Chemistry and Courant Institute of Mathematical Sciences
New York University, New York NY 10003
Glenn J. Martyna
Dept. of Chemistry,
Indiana University, Bloomington, IN 47405
One of the ``grand challenges'' is to develop methodology capable of providing an accurate description of conformational equilibria in systems with rough energy landscapes via computer simulation. If met, many important problems, most notably protein folding, could be impacted. Here, a significant step forward is made by combining with molecular dynamics a novel variable transformation designed to warp configuration space such that barriers are reduced and basins are enlarged while equilibrium properties are preserved. Thus, barriers are readily crossed and sampling is very efficient.
One of the premier challenges in molecular simulation is to determine accurately and efficiently conformational equilibria of systems described by complex potential surfaces that possess high energy barriers separating important minima or basins of attraction. Indeed, rough energy landscapes are ubiquitous in physical, chemical and biological problems, including protein folding. Despite recent advances [, ] this important class of problems cannot, yet, be adequately addressed by molecular simulation.
Current approaches include umbrella sampling, guiding potentials, multicanonical algorithms, and rate-enhancing schemes. The first two techniques reweight the phase space distribution in order to accelerate sampling at the price of requiring an a posteriori correction. Thus, the efficiency scales exponentially poorly with the number of degrees of freedom involved in the reweighting factor. Parallel tempering methods introduce a series of M simulations at different temperatures which are permitted to interchange. Thus, the cost of the calculation is M-fold larger, and the number of temperature interchanges accepted scales exponentially poorly with the system size. Other methods that rely on harmonic transition state approximations do not, in general, preserve equilibrium properties.
In this letter, a novel approach to the conformational sampling
problem that is free of the above difficulties is presented.
The method is based on an exact reformulation
of the classical statistical mechanical
partition function via a nonlinear transformation of
the integration variables. The transformation is constructed so
as to effect a warping of the configuration space
wherein barrier regions are shrunk, barrier heights
are lowered, and attractive basins stretched.
The result is a smoother surface, which allows barriers to
be easily traversed without i) altering equilibrium properties of the system,
ii) requiring an a posteriori reweighting,
iii) relying on harmonic approximations,
or iv) scaling exponentially poorly with the number
of transformed degrees of freedom.
The technique is combined, here, with molecular dynamics (MD)
to form a novel conformational sampling scheme, although
it could also be be used with Monte Carlo.
Moreover, the new approach could be combined with complementary
transformations [] that help overcome entropic barriers as well
as employed to enhance the efficiency of other techniques such as
parallel tempering.
Applications of the new method to model and complex systems as
large as 400-mer alkane chains are considered in order to demonstrate
its favorable scaling with system
size and its enhancement in sampling efficiency
(a factor of 
over standard MD for the long chains, as will be seen shortly).
Consider, first, a classical particle with momentum, p, mass, m, and coordinate, x, in a one-dimensional potential, V(x), which possesses two stable minima, separated by a large energy barrier. The canonical partition function is
where h is Planck's constant, and
. The spatial probability distribution function,
,
appearing in Eq. (1), could be sampled using a classical
Hamiltonian, 
, in conjunction with
thermostatted MD [].
However, if the potential
barrier is too high, barrier crossing events will be rare, and extremely
long trajectories will be needed to generate an adequate sampling of
the configuration space.
Since the variable, x, is just an integration variable, a transformation
to a new variable, u=f(x), can be performed arbitrarily without
altering the partition function in Eq. (1), assuming that the
inverse 
exists and is single-valued. Upon changing variables
and exponentiating the Jacobian factor g'(u),
the partition function becomes
The key feature of Eq. (2) is the appearance of an effective potential
that can be adjusted via the transformation to be smoother than V(x). It is important to note that this approach completely eliminates the need for a posteriori reweighting factors required by the umbrella and guiding potential methods and does not introduce harmonic approximations. It is, therefore, not equivalent to any other existing schemes.
A spatial warping transformation is achieved via the following change of variables:
Here, 
is an arbitrary reference potential,
is an arbitrary point, and c is
a constant. Since the integrand in Eq. (4)
is nonnegative and u
is a monotonically increasing function of x, a unique
inverse x=g(u) exists. The Jacobian of the transformation is
Substituting Eq. (4) and the Jacobian into Eq.\
(2)
gives an effective potential
.
If is taken to be equal to the bare
potential, V(x), in the barrier region and zero outside
and is taken to be the left minimum ( 
),
then 
is barrier free (see Fig. 1(a)). Thus, the
Jacobian weights the ``u'' space so as to reduce the barrier
region to a negligibly small volume without altering the partition function.
We refer to the new approach as
REPSWA (REference Potential Spatial Warping Algorithm).
Figure 1: (a) Schematic of REPSWA for a double-well potential.
(b,d) x(t) in the double well potential
( 
) without(b)/with(d) REPSWA.
(c,e) Probability distribution without(c)/with(e) REPSWA (dashed
lines) and analytical result (solid line).
The reformulated partition function, Eq. (2), can be sampled by a canonical MD approach based on a Hamiltonian,
forming the REPSWA-MD method.
Since the transformation is not canonical, 
generates trajectories
which sample phase space more effectively than those of the
original Hamiltonian,
That is, the true dynamics, in which barrier crossing is, by definition, a rare event, is not, nor is intended to be, preserved.
Since the integral in Eq. (4) cannot generally
be evaluated analytically,
the function, 
, is expressed a finite
basis of integrable functions, thereby defining a transformation 
,
so that the integrals can be performed.
The function, 
,
will be very close to f(x) in Eq. (4)
for a large enough basis set.
This is a key step that allows the REPSWA-MD to be applied beyond the harmonic
approximation.
Clearly, the transformation, also exactly preserves
the partition function.
The REPSWA-MD technique is first applied to a one-dimensional analytically solvable system in order to demonstrate its validity and performance. A double well potential of the form
is selected, and canonical MD simulations are performed
for a high barrier, .
The equations of motion for 
are coupled to a
Nosé-Hoover chain thermostat [] to effect the canonical sampling.
Runs of length of 10 
steps
were performed using a time step of 
( 
,m=1 and a=1).
Figure 1 shows trajectories x(t)=g(u(t)) with and without
the REPSWA transformation as well as the
probability distribution function, P(x).
The ability of the REPSWA technique to
reproduce the correct distribution, P(x),
demonstrates that the partition function is perfectly preserved.
In contrast,
ordinary canonical MD based on H
is unable to generate
the correct distribution within steps indicating, quantitatively,
an increase in efficiency of over 6 orders of magnitude using REPSWA.
Next, the nontrivial extension of the REPSWA method to enhance sampling of
multidimensional potential surfaces is discussed. In particular,
chain molecules possess many dihedral barriers
in the range of 5-10 
, which, as the previous
example has demonstrated, are traversed infrequently, hindering
conformational sampling. It will first be shown
how to treat the commonly used ``united atom'' (UA)
molecular models (e.g. CH 
-like moieties are merged into
pseudoatoms) [], and
then the additional steps needed
to treat all-atom models [] will be discussed.
Dihedral barriers can be removed in united atom models
by applying the REPSWA transformation
(cf. Eq. (4)) directly to the dihedral angle. However, since
dihedral angles are not explicit coordinates
in MD, additional steps are needed in order
to apply the technique. Consider a UA representation of butane,
and assume the connectivity is C 
-C 
-C 
-C 
, where
C 
, the ith pseudoatom, has Cartesian position 
.
The transformation scheme is:
(i) the vector 
is rotated/translated into a frame in
which 
is at the origin and 
lies along the
z axis, (ii) the new 
is then resolved into spherical polar
coordinates, 
, (iii) REPSWA is applied to the
azimuthal angle to generate
, (iv) a Cartesian 
is then created using
, (v) is placed back into the original frame
by inverting the transformations (i)
and (ii) to obtain a transformed set of Cartesian coordinates.
The constant c is adjusted to ensure that 
.
The reference potential in step (iii) is chosen to be equal to the
dihedral potential
between the two gauche conformations and zero otherwise.
The forces on the 
are obtained by the chain rule. This procedure
allows for facile transformation of all the dihedral angles
in an N-pseudoatom chain by moving sequentially down the chain.
This procedure forms an upper triangular Jacobian matrix whose determinant is the
product of the
exponentials of the reference potentials, which cancel the desired barriers,
while preserving the partition function.
Both the
coordinate and force transformations can be carried out in 
operations.
The behavior of the dihedral REPSWA scheme is studied
using the Ryckaert-Bellemans model []
for a 400-mer chain at T=300K (ignoring intermolecular
Lennard-Jones interactions to avoid freezing the chain).
It should be noted that guiding potentials/umbrella sampling
could not be used to enhance the sampling of all 397 dihedral angles.
The molecule was simulated using canonical molecular dynamics
in the gas phase with a time step of 0.5 fs for 10 steps
with and without the REPSWA dihedral angle transformation.
Figure 2 demonstrates the large improvement in sampling
efficiency given by REPSWA.
The evolution of the end-to-end distance, 
( 
Å, N=400)
(Fig. 2(c), (f)) indicates that
the fluctuations in this quantity are virtually nonexistent
without REPSWA.
Quantitatively, the number of dihedral angle flips per torsion
increases by a factor of 30, and the correlation time for the
relaxation of 
decreases by 6 orders of magnitude!
Note, the REPSWA calculation was found to be in good agreement with the
prediction ( 
)
of a three state rotational-isomeric state (RIS) model [],
parameterized using our potential.
Figure 2: (a,d)
Dihedral angle,
,
in the UA 400-mer without(a)/with(d) REPSWA.
(b,e) Number of dihedrals, 
,
that have crossed barriers 
times without(b)/with(e) REPSWA.
(c,f) End-to-end distance ( , 
Å)
without(c)/with(f) REPSWA.
Last, it is shown how REPSWA can be used to simulate
``all-atom'' models of chain molecules with high efficiency
(see Fig. 3(a)).
In all-atom models, the backbone C -C -C -C
dihedral angle is strongly coupled to the
H -C -C -C , H -C -C -C ,
and H -C -C -C dihedrals due to
the 
hybridization of the carbon atoms (represented
by harmonic bond and bend angle potentials).
The rigidity caused by the hybridization requires groups of atoms
to be transformed as a rigid units. Groups of three atoms
(e.g. C , H and H 
form such a group in Fig. 3(a)) are
selected, and a ``primary'' dihedral angle (e.g. H -C -C -C ), defined.
A set of ``secondary'' dihedrals is formed by collecting
dihedral angles involving atoms both in and to the left of the group
(e.g. eight). The reference potential is then taken to be
where 
and
are the primary and secondary dihedral potentials, respectively,
is a function that switches off
the potential outside the region between the two gauche conformations,
and 
is a set of constants chosen based on the
ideal differences between the secondary and primary dihedral angles
for geometry. In a simulation, the actual secondary
dihedrals will generally not differ much from 
.
The full REPSWA transformation procedure
is then the same as that for the
united atom case, except that three vectors 
,
and 
are
rotated about the C -C axis.
The performance of the group REPSWA transformation scheme is tested on
an all-atom 400-mer using the CHARMM22
force field [].
The simulation details are the same as in the united
atom case (see above).
Figures 3(b) and 3(c)
show the histogram of backbone dihedral angle ``flips'', again,
without and with REPSWA, respectively; the insets
show the evolution of a typical backbone dihedral for the two cases.
It can be seen
that the use of REPSWA leads to a much greater dihedral flipping rate
by a factor of 30, and, hence, an improved sampling of
conformational space by a factor of based on
(Figs. 3(d) and 3(e)).
Again, the REPSWA calculation was found to be in good agreement with the RIS
model [] prediction ( 
).
Figure 3: (a) Schematic of group REPSWA for C H 
.
(b,d) Number of dihedrals, ,
that have crossed barriers times without(b)/with(d) REPSWA;
inset: .
(c,e) End-to-end distance ( , Å)
without(c)/with(e) REPSWA.
In conclusion, a novel method for enhancing conformational sampling simulations of complex molecular systems, REPSWA-MD, has been introduced. REPSWA-MD smooths a rough potential energy surface via a nonlinear variable transformation and, thus, does not alter the equilibrium properties of a system. Unlike guiding potentials or umbrella sampling, the method does not require a posteriori reweighting and can be used to enhance sampling in a very large numbers of degrees of freedom without exponentially poor scaling. A dramatic improvement in the rate of barrier crossing events and the number of statistically independent conformations sampled has been demonstrated (e.g. 6 orders of magnitude for a 400-mer chain). The new technique can easily be combined with complementary transformations [] that help overcome entropic barriers as well as with parallel tempering schemes. Our current efforts are focused on applying REPSWA-MD to proteins and other large biomolecules.
This research was supported by PRF 32139-AC and NSF CHE-96-5015 (G.J.M.) and by Research Corp. RI0218, NSF CHE-98-75824, and an NYU Whitehead Fellowship in Biomedical and Biological Sciences (M.E.T.) and NSF-EIA-0081307 (G.J.M and M.E.T.).
