J. Phys. Chem. B 2000, 104, 6877-6883
6877
A Stochastic Dynamics Simulation Study Associated with Hydration Force and Friction
Memory Effect
Ben Zhuo Lu,†,‡ Cun Xin Wang,*,† Wei Zu Chen,† Shun Zhou Wan,‡,§ and Yun Yu Shi‡
Center for Biomedical Engineering, Beijing Polytechnic UniVersity,
Beijing 100022, People’s Republic of China, and School of Life Science, UniVersity of Science and
Technology of China, Hefei 230026, People’s Republic of China
ReceiVed: December 10, 1999; In Final Form: March 22, 2000
A new simulation approach for combining hydration force with generalized Langevin dynamics is developed
in this paper. The exponential model is taken for the friction kernel. The hydration force determined by the
boundary elementary method is taken into account as the mean force terms of the solvent, including the
Coulombic interactions with the induced surface charge and the surface pressure of solvent. All simulations
were performed on cyclic undecapeptide cyclosporin A (CPA). The simulation results obtained using the
new method were analyzed and compared with those obtained using other methods, such as molecular dynamics
simulations, generalized Langevin dynamics simulations, and conventional stochastic dynamics simulations.
We found that the results obtained with the new method presented in this study show obvious improvements
over the other simulation techniques and that the hydration force and friction relaxation together contribute
to this improvement.
Introduction
Molecular dynamics (MD) simulations have been widely used
to study the structural and dynamic properties of molecular
systems.1,2 However, there exist at least two limitations on this
approach:1 the approximation in the potential energy functions
and the lengths of the simulations. The first limitation leads to
systematic errors, and the second to statistical errors. Some
previous studies have been performed to improve the form of
potential functions.3-5 To prolong the simulation time, the use
of stochastic dynamics (SD) based on the Langevin equation is
a recommended approach in which only the relevant portion of
the molecule is considered explicitly and the remainder of the
system, such as the solvent, serves to provide an effective
potential, a friction drag, and a heat bath.7 The friction memory
effect is taken into account in generalized Langevin dynamics
(GLD), which is a more proper approach when the solvent and
solute particles have similar sizes and masses.
A potential of mean force of solvation that describes the
average solvent effect on the solute degrees of freedom has
important influence on the solute conformation.5 However,
conventional SD and GLD simulations have usually omitted
the extra mean force terms of the solvent such as the Coulombic
interaction with the induced surface charge and the boundary
pressure exerted by the solvent.7,8,10,12 Many efforts have been
made to incorporate solvent effects into molecular mechanics
or MD simulations using the Poisson-Boltzmann equation.8-11
Recently, Wang and co-workers succeeded in incorporating the
hydration force determined by the boundary elementary method
* Author to whom correspondence should be addressed. E-mail: cxwang@
bjpu.edu.cn.
† Beijing Polytechnic University.
‡ University of Science and Technology of China.
§ Present address: Laboratoire de Chimie Biophysique, Universite Louis
Pasteur ISIS 4, rue Blaise Pascal, 67000 Strasbourg, France.
(BEM)13 into conventional SD simulations (SDBEM).4 Then,
Wan and co-workers used an exponential model for the friction
kernel and developed the leapfrog algorithm for numerical
integration of generalized Langevin dynamics.12 Until now,
however, no one has reported an efficient approach to incorporating both hydration force and friction memory into SD
simulations for biomolecules. The present work is a continuation
of the work developed by Wang et al.4 The goal of this work is
to develop a new simulation approach to GLD with a potential
of mean force calculated by the BEM (GLDBEM) and to clarify
how the hydration force and friction memory affect the
conformational and dynamical behavior of the system in the
GLDBEM simulation. Because we found that the original GLD
algorithm by Wan et al.12 may be subject to mathematical
difficulty in actual calculations dealing with the stochastic terms
in the simulation, the related parts of the GLD programs have
been modified in the present study to solve such problems. To
compare our simulation data with those obtained in the previous
work, we have selected the same system, cyclosporin A (CPA),
as our simulated system for examining the rationality and
reliability of our GLDBEM approach.
Theory
For a system of N atoms, the GLD technique is based on the
following generalized Langevin equation (GLE):
miẍi(t) ) Fi({xi(t)}) - mi
∫0t γi1(t - τ) x̆i(t) dτ + Ri1(t)
(1)
where xi is the coordinate component and mi the mass of the
ith atom, Fi({xi(t)}) represents the system force, γi1(t - τ) is
the friction kernel, and Ri1(t) denotes the random force, which
satisfies the fluctuation dissipation theorem14
10.1021/jp994349o CCC: $19.00 © 2000 American Chemical Society
Published on Web 06/28/2000
Lu et al.
6878 J. Phys. Chem. B, Vol. 104, No. 29, 2000
〈Ri1(0)Ri1(t)〉 ) mikBTγi1(t)
(2)
where the angular brackets indicate an ensemble average, kB is
the Boltzmann constant, and T is the absolute temperature.
The Laplace transform of the higher-order memory function
γin(t) can be expressed by the continued fraction equation15
γ̃1(s) )
γ1(0)
s + γ̃2(s)
γ1(0)
)
(3)
γ2(0)
s+
s+
γ1(0)
s + 1/τ0
(4)
γi1(t) ) γi0 exp(-t/τ0)
(5)
where τ0 is the relaxation time of the friction kernel and γi0 is
its initial value. Therefore, the exponential model for the friction
kernel in eq 5 is also taken in this work as a simple approximate
form. In this case, the leapfrog algorithm16 can be used for
numerical integration of the GLE with the exponential model
for the friction kernel, and the numerical integration formula
and procedure of GLD can be found in ref 12.
Here, we should point out that some modifications are made
in our present work. The leapfrog algorithm was adopted in
our work for the integration of eq 1. Therefore, there are two
sets of random variables of Wn(-∆t/2), Xn-1/2(∆t/2), Vn(-∆t/
2) and Wn(∆t/2), Xn+1/2(-∆t/2), Vn(∆t/2) in the leapfrog
algorithm, in which W is a random variable with a Gaussian
distribution and X and V are random variables appearing in the
integration expressions for the position and velocity, respectively.12 Each set of random variables obeys a trivariate Gaussian
distribution17
1
(2π)3/2|B|
mγ0kBT
×
τ0
(ξτ0)2 ∆t/2 - 2τ0(1 - exp(-∆t/2τ0)) +
]
τ0
(1 - exp(-∆t/τ0)) - 2ξτ0[∆t/2 2
τ0(1 - exp(-∆t/2τ0)) - ξ-1(1 - exp(-ξ∆t/2)) +
(ξ + 1/τ0)-1 × (1 - exp(-∆t/2τ0) exp(-ξ∆t/2))] +
}
where
ξ = γ0τ0[1 - exp(-∆t/t0)]
}
1
exp - (S - M)T B-1 (S - M) (6)
2
where S is the random vector whose components are the random
variables of s1, s2, and s3 and M is the corresponding mean
vector, M ) 〈S〉. B is the covariance matrix [σij], σij ) 〈sisj〉1/2,
and |B| its determinant. In practical calculations, we have found
that the determinant |B| of the positive determined matrix B
may be given by a negative value because of the calculating
precision of the simulation. To solve this problem, we used
power expansions for the matrix elements. For instance, the
matrix element of B for the first set of random variables, σ22,
has the form (see eq B2 in Appendix B of ref 12)
(8)
When we used the power expansions for the matrix elements,
the power expansion form of σ22 can be given by
2
(∆t/2)〉 )
σ22 ) 〈Xn-1/2
[
(
mξ2 1 -
1
ξτ0
)]
-2
(
[ (
(
)
mγ0kBT (∆t)5 4
ξ3 ξ2
ξ -2 + 2 τ0
640
τ0 τ
)
0
(∆t)6 5 ξ4 ξ3
ξ2
ξ - - 2+ 3 +
2304
τ0 τ
τ
0
0
)]
(∆t)7
ξ3
ξ3
ξ2
ξ5
5ξ6 - 3 - 4 2 - 3 2 + 5 4
64512
τ0
τ
τ
τ
0
0
0
+ ... (9)
A similar operation on the other elements of matrix B can be
done in the same way. When we used this procedure in this
work, the problem with calculating the value of the determinant
|B| mentioned above was solved.
The system force Fi({xi(t)}) in eq 1 can be written as
Fi(t) ) - ∂Vmean({xi})/∂xi
(10)
where Vmean is the system potential of mean force. In the present
work, Vmean can be expressed as
Vmean({xi}) ) Vint({xi}) + Vsol({xi})
(11)
where Vint refers to the internal potential of solute atoms in the
system and Vsol is the potential of mean force from the solvent.
Hence, the system force can be obtained using
Fi ) -
×
1/2
{
-2
(7)
Thus, the second-order of the memory function is γ2(t) ) 1/τ0
δ(t), and the original function of the Laplace transform γ̃1(s)
has an exponential form
f(s1,s2,s3) )
{ [
)]
1
ξτ0
∆t/2 - 2ξ-1[1 - exp(-ξ∆t/2)] + (2ξ)-1[1 - exp(-ξ∆t)]
γ3(0)
s + ...
where γn(t) is the nth order memory function and γ̃n(s) the
Laplace transform of γn(t). If the continued fraction in eq 3 is
truncated at n ) 2, the Laplace transform of γ1(t) can be written
as
γ̃1(s) )
[ (
2
(∆t/2)〉 ) mξ2 1 σ22 ) 〈Xn-1/2
∂Vint({xi}) ∂Vsol({xi})
sol
) Fint
i + Fi
∂xi
∂xi
(12)
where Fint
i is the force due to the simulated atoms of the system
and Fsol
i is the hydration force due to the environment solvent.
In this work, the hydration force Fsol
i is calculated using the
classical continuum model of electrostatic interactions, which
includes two parts: one is the Coulombic interaction with the
induced surface charge, and the other is a purely mechanical
boundary pressure of the solvent.8,10 The hydration force is
calculated by the BEM in this study. Because the detailed
procedure of the BEM can be found in many other works,13,18-20
here, we only give the main equations used in this study.
The interior potential φi and the exterior potential φe are the
solutions of the Poisson equation and the linearized Poisson
Boltzmann equation, respectively.
Hydration Force and Friction Memory Effect
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6879
1
∑qk δ(rp - rk)
D k
∇2φip ) -
(13)
charged atom and the uncharged atoms on the surface can be
calculated by the expression8,10
[ (
i
∇2φep
) k2φe(rp)
(14)
where Di is the dielectric constant of the molecular interior, rp
is a point inside or outside the molecule, qk is the kth point
charge at rk, and k in eq 14 is the Debye inverse screening
length. The solutions of eqs 13 and 14 can be written in the
integral forms by using Green’s second theorem on eqs 13 and
14.
φip
)
[
∂φik
I Gpk ∂n - φik
S
φep )
∂n
[
1
dAk +
]
∑k qkGpk
Di
∂φek
∂upk
+ φek
-upk
dAk
∂n
∂n
I
S
]
∂Gpk
(15)
(16)
where the integration is taken over the entire molecular surface,
k is a point on the surface, p is a point inside or outside the
molecule, and n is the outward unit normal to the surface. The
fundamental solution of eqs 13 and 14 can be expressed as
Gpq )
1
4πrpq
(17)
upq ) exp(-krpq)/4πrpq
(18)
where rpq is the distance between points p and q. When the
jump discontinuity of the double-layer potential at the boundary
S is considered,18 eqs 15 and 16 can be written as
1
2
φip )
I
S
[
Gpk
1 e
φ )
2 p
∂φik
∂n
[
I
S
- φik
-upk
]
∂Gpk
∂n
dAk +
1
Di
]
∑k qkGpk
∂φek
∂upk
dAk
+ φek
∂n
∂n
(19)
(20)
The surface potential and its normal derivatives can be evaluated
from eqs 19 and 20 plus the following two boundary conditions:
φi ) φe
(21)
i
Di
e
∂φ
∂φ
) De
∂n
∂n
(22)
where De is the dielectric constant of the solvent.
The electrostatic component of the hydration force acting on
charge qi is calculated with the expression
Fpol
i (R) )
I
S
qiσ(r)(R - r)
|R - r|
3
d 2r
(23)
where the integration is carried out over the entire molecular
surface S and σ(r) is the density of the induced polarization
charge at an arbitrary point r on the molecular surface S. The
relation between the σ(r) and the exterior electric field Ee(r) is
given by19
σ(r) )
(
)
1 Di - D e
Ee(r)‚n(r)
4π
Di
(24)
The purely mechanical boundary pressure acting on both the
p(ri) )
) | |]
De - D i
1
+ E0
(D - Di) (Ee‚n)2
8π e
Di
2
(25)
Thus, the force of solvent pressure on the ith atom is given by
)
Fpress
i
∫S - np(r) d2r
(26)
where the integration is performed only over the portion of
solvent accessible surface associated with the ith atom.
Computational Procedures
The hydration force and friction memory effect were incorporated into the conventional SD simulation technique
simultaneously in this work by linking the BEM programs
(MACBEM)13 with the GLD program12 based on the GROMOS
package.21 The program for calculating the hydration force on
a charge and the solvent pressure based on the BEM were joined
with the subroutine program FORCE in the GROMOS package.
In the present work, the triangulation procedure in the BEM
was similar to the method developed by Juffer and co-workers.20
The molecular surface was defined as the center of the probe
rolling around the molecule. The probe radius is 0.16 nm. The
number of total triangles for CPA was 320.
In GLD, the atomic friction kernels satisfy the relation
∫0∞ γ(t) dt ) γτω
(27)
where γτ is the friction constant, which was assigned a value
of 91 ps-1 to represent the total solvent effect of water
molecules,7 and ω is an atomic accessible area weight factor,
which is due to the degree of solute-solvent interaction. The
calculation of ω usually uses the approximate expression
ω ) max(0, 1 - Nnb/Nnbref)
(28)
where Nnb is the number of neigbour atoms of the concerned
atom within a sphere of radius Rnbref and Nnbref is a reference
number. In this work, we chose Nnbref ) 10. The parameter Rnbref
was set to 0.3 nm. The relaxation time, τ ) 0.1 ps, was used
for the friction kernel,12 and the initial friction coefficients, γi0,
were taken to be 910 ps-2, weighted with the accessible area
factor.
The GROMOS force field21 was used for all simulations. The
initial structure of the CPA molecule used was the X-ray
structure,22 which is shown schematically in Figure 1. Nonpolar
hydrogen atoms were included in the carbon atoms (united atom
approach), while polar hydrogen atoms were treated explicitly.
The system contained 90 atoms in total. The time step was taken
to be 2 fs for integrating the equations of motion. A temperature
bath and a pressure bath were applied to keep the system at
300 K and 1 atm.24 All bond lengths were kept rigid using the
SHAKE algorithm with a tolerance of 10-4.25
The length of the simulation is one of the main concerns with
protein dynamic simulations. Some recent studies suggested that
a simulation length of 500 ps is acceptable. Daggett and Levitt
pointed out that any solution simulation under 50 ps in duration
is probably not sufficiently equilibrated for one to draw any
conclusions about the behavior of proteins.23 Because the
previous MD, SD, and GLD simulations4,6,7,12 were only
performed for 40-60 ps, they seemed to be not long enough.
In our work, after energy minimization, a 100-ps conventional
Lu et al.
6880 J. Phys. Chem. B, Vol. 104, No. 29, 2000
TABLE 1: RMS Atomic Positional Fluctuationsa and Mean
Anisotropy for Atomic Motion in CPA
atom type
Figure 1. Schematic structure of cyclosporin A. The residue name of
MeBmt is abbreviation for 4-(E-2-butenyl)-4,n-dimethylthreonine. The
broken lines show the hydrogen-bonding network with high occupancies. A few of the atoms that contribute to the formation of the hydrogen
bonds appearing in Table 2 are labeled, such that H refers to the
hydrogen atom bonded to the nitrogen atom on the backbone and OH
to the hydrogen atom in the hydroxyl group in 1MeBmt.
MD
all atoms
CR atoms
MeLeu Cβ
MeLeu Cγ
MeLeu Cδ
0.084
0.052
0.067
0.099
0.141
all atoms
CR atoms
0.38
0.43
a
SD
Results and Discussion
Position Fluctuation and Deviation. The root-mean-square
(RMS) positional fluctuation and the RMS deviation of the
generated structure from a reference structure as a function of
time are important parameters for examining the molecular
flexibility and the convergence of the simulations. Figure 2
shows the time-dependent RMS deviation of atomic positions
in the GLDBEM simulation (after the 100-ps equilibration)
relative to the X-ray crystal structure. As shown in Figure 2,
the RMS shift in position has a stable value of 0.08-0.21 nm
for all atoms and 0.04-0.14 nm for CR atoms during the 500ps simulation, which corresponds with the previous simulation
work.23 It is found that there are two small periodical peaks in
Figure 2. The slight instability of the positional fluctuation can
also be found in another study,23 in which various patterns of
positional fluctuation and deviation were reported and it was
mentioned that the deviation could result from many different
factors. Because the system can occur periodically with a slight
reasonable deviation, the system seems to be in an equilibrium
state after 100 ps of the GLDBEM simulation.
Table 1 lists the RMS atomic positional fluctuations and the
mean anisotropy of atomic motion in CPA for the MD, SD,
0.35
0.39
GLDBEM
0.078
0.049
0.076
0.102
0.138
0.36
0.37
Fluctuations are in nanometers.
TABLE 2: Frequencies of Intramolecular Hydrogen Bondsa
Obtained from Different Simulation Techniques in CPA
donor
a
SD simulation was performed for system equilibration. Then,
a 500-ps GLDBEM simulation was carried out, and the
trajectories were saved every 25 time steps. For comparison,
we also performed MD, GLD, and SD simulations again for
500 ps each. The trajectories from the last 400 ps of all
simulations were used for analysis.
SDBEM
Anisotropy
0.38
0.37
0.45
0.38
acceptor
1MeBmt O-H 10Meleu
2Abu
N-H 5Val
N-H 11MeVal
5Val
N-H 2Abu
N-H 3Sar
7Ala
N-H 5Val
N-H 11MeVal
8Ala
N-H 6MeLeu
Figure 2. Root-mean-square deviation of atomic positions relative to
the X-ray structure versus GLDBEM simulation time.
GLD
Positional Fluctuations
0.065
0.074
0.075
0.040
0.044
0.044
0.061
0.068
0.075
0.081
0.094
0.104
0.122
0.132
0.137
MD SD GLD SDBEM GLDBEM
O 0
O 16
O 4
O 67
O 2
O 0
O 13
O 8
2
32
4
76
1
1
51
44
3
36
4
67
4
2
52
43
2
26
7
64
2
3
50
37
1
27
4
52
2
3
45
36
Frequencies are in percentages.
GLD, SDBEM, and GLDBEM simulations. The RMS positional
fluctuation of our GLDBEM simulation is 0.078 nm for all
atoms, which is larger than those of the SD (0.065 nm), GLD
(0.074 nm), and SDBEM (0.075 nm) simulations and closer to
that of the MD simulation (0.084 nm), which gives the largest
positional fluctuation. This result is consistent with the recent
simulation work by Fraternali and van Gunsteren.5 Note that
the value of the RMS fluctuations for the longer MD simulation
(500-ps) is larger than that obtained from the previous shorttime (50-ps) MD simulation (0.052 nm for all atoms and 0.036
nm for CR atoms).6 Daggett and Levitt have reported similar
results obtained from different lengths of MD simulations for
the different systems.23 Generally, the RMS positional fluctuation in long simulations is larger. In the present work, the RMS
positional fluctuation in the 500-ps MD simulation is the largest.
This means that the CPA molecule shows more flexibility in
the MD simulation in solution than in the other simulations
without explicit solvent. As shown in Table 1, the values of
the RMS fluctuations of the GLDBEM simulation are larger
than those of the GLD simulation, and the RMS fluctuations of
the SDBEM simulation are larger than those of the SD
simulation. This means that the GLDBEM and SDBEM
simulations can partially reflect the influence of the hydration
force determined by the BEM. Comparing GLDBEM with
SDBEM, and GLD with SD, it is found that the RMS
fluctuations of the former simulations are larger than those of
the latter simulations. This reflects the fact that the friction
memory effect imposes an obvious influence on the motion of
the molecular atoms. The mean anisotropy in the atomic motion
in Table 1 is defined as the ratio of the shortest to the longest
principal axis of the anisotropic fluctuation ellipsoids. The values
for anisotropy are comparable for all five simulations.
Hydrogen-Bonding Analysis. Table 2 reports the results of
a hydrogen-bonding (H-bonding) analysis of CPA obtained from
MD, SD, GLD, SDBEM, and GLDBEM trajectory data for the
last 400 ps. The criteria used to determine a H-bond are purely
geometric: for each coordinate set, every potential donoracceptor pair is tested and considered to form a H-bond if the
hydrogen-to-acceptor distance is less than 0.25 nm and the
Hydration Force and Friction Memory Effect
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6881
donor-acceptor angle is larger than 140°. The frequencies of
H-bonds are determined from the occurrences registered on the
simulation trajectory frames.
In Table 2, eight H-bond donor-acceptor pairs and their
occupancies are listed for different simulation approaches. The
H-bond occupancies among the pairs 1MeBmt-10MeLeu,
2Abu-11MeVal, 5Val-3Sar, and 7Ala-5Val obtained from
the GLDBEM simulation are close or equal to the results from
MD simulation. The H-bonds between the pairs 7Ala-11MeVal
and 8Ala-6MeLeu in the four stochastic simulations have much
higher occupancies than those in the MD simulation. It is found
that the residues of 7Ala, 11MeVal, and 8Ala are located at the
edge of the loop region. Because the structure of the loop region
of CPA in aqueous solution is more open than that of the
β-pleated sheet, the internal H-bonds between the pairs 7Ala11MeVal and 8Ala-6MeLeu in the MD simulation have lower
occupancies because of the competition of water molecules for
hydrogen-bond donors and acceptors in the solute.
As shown in Table 2, for the H-bonding patterns among the
pairs 1MeBmt-10MeLeu, 2Abu-5Val, 5Val-2Abu, 5Val3Sar, 7Ala-11MeVal, and 8Ala-6MeLeu, the results obtained
from the GLDBEM simulation show lower occupancies compared with the data from GLD simulation. In Table 2, it is also
found that the H-bond occupancies of the SDBEM simulation
are lower compared with those of the SD simulation for the
H-bonds between the pairs 2Abu-5Val, 5Val-2Abu, 7Ala11MeVal, and 8Ala-6MeLeu. This means that the GLDBEM
and SDBEM simulations can partially reflect the mean force
of the solvent, which imposes an influence on the conformation
of the H-bond.
In addition, there are four pairs of H-bonds (1MeBmt10MeLeu, 2Abu-11MeVal, 7Ala-11MeVal, and 8Ala-6MeLeu) in Table 2 in which the occupancies of the H-bonds
obtained from the GLDBEM simulation are closer to those from
the MD simulation than to those from the SDBEM simulation.
This indicates that the friction memory effect in the GLDBEM
simulation also has an influence on the H-bonding network.
Dynamical Behavior. The friction and stochastic terms in
the GLE and the hydration force both influence the dynamical
properties that can be reflected by the time correlation functions.
The time evolution of the atomic positional autocorrelation
function has been examined using the following autocorrelation
function:
c(t) )
〈∆r(t) ∆r(t + τ)〉τ
〈∆r 〉
2
)
〈∆r(0) ∆r(t)〉
〈∆r2〉
(29)
where the brackets 〈‚‚‚〉 represent an average over the simulation
time, ∆r(t) ) r(t) - 〈r〉 , and 〈r2〉 ) 〈(r - 〈r〉)2〉.
In this study, we have selected two backbone atoms and two
side-chain atoms of CPA for analyzing the positional autocorrelation functions for different simulation approaches. Figure
3a-d shows the autocorrelation functions for those four atoms
of CPA, in which the backbone atom 7Ala-CR and the sidechain atom 4Meleu-Cδ1 are in the β-pleated sheet and the
backbone atom 10Meleu-CR and the side-chain atom 10MeleuCδ1 are in the loop region. As shown in Figure 3a and b, the
correlation functions for backbone atoms 7Ala-CR and 10MeleuCR have a rapid initial decay followed by a linear decrease for
all simulations. However, the correlation functions of the GLD
and SD simulations in Figure 3b decay even more quickly than
those of the other simulations. The correlation function obtained
from the GLDBEM simulation in Figure 3b is closer to that
obtained from the MD simulation compared with other simula-
Figure 3. Atomic positional fluctuation autocorrelation functions for
four atoms of CPA: (a) 7Ala-CR, (b) 10MeLeu-CR, (c) 4MeLeu-Cδ1,
and (d) 10MeLeu-Cδ1. The bold solid lines show the MD simulation,
the thin solid lines show the GLDBEM simulation, the dashed lines
show the SDBEM simulation, the dash-dotted lines show the GLD
simulation, and the dotted lines show the SD simulation.
Lu et al.
6882 J. Phys. Chem. B, Vol. 104, No. 29, 2000
Figure 4. Atomic positional fluctuation autocorrelation functions for
10MeLeu-Cδ1 of CPA for the GLDBEM simulation obtained by
averaging for five different simulation times (20, 40, 100, 200, and
400 ps).
tions, and the same case can be found in Figure 3a for atom
7Ala-CR in the early period of 10 ps. This indicates that our
GLDBEM simulation for these two backbone atoms gives
meaningful results for the correlation functions compared with
those of other stochastic simulations. For the atom 10MeLeuCδ1 (Figure 3d), the correlation functions obtained from the
GLDBEM and SDBEM simulations are very consistent with
those obtained from the MD simulation in the early 10-ps period.
As shown in Figure 3c for the atom 4MeLeu-Cδ1, the correlation
functions of the SDBEM and GLDBEM simulations show even
slower decay than those for the MD and other simulations. From
the structure analysis mentioned above, the large RMS positional
fluctuations of the atom 4MeLeu-Cδ1 can be found in the
GLDBEM (0.178 nm) and SDBEM (0.214 nm) simulations.
This could be the main reason for the slow decay of the
correlation functions shown in Figure 3c for the SDBEM and
GLDBEM simulations.
It is interesting to find that the correlation functions of the
GLDBEM simulations in Figure 3a-c are closer to those of
the MD simulations than to those of the SDBEM simulation.
This reflects the influence of the friction memory effect on the
atomic dynamical properties of the system. In addition, the
correlation functions in Figure 3a-d for the GLDBEM and
SDBEM simulations show generally longer relaxation times than
those for the SD and GLD simulations. This indicates that the
mean force of hydration in the GLDBEM and SDBEM
simulations also impose somewhat of an effect on the relaxation
time.
It is worth noting that, in our present work, the correlation
functions were obtained by averaging on a scale of a 400-ps
simulation time, whereas in the previous work,4,6,7,12 the
averaged simulation period was only 40 ps. To understand the
difference in the correlation functions obtained from the different
simulation time scales, we have analyzed the atomic positional
correlation function for the GLDBEM simulation using different
simulation times. The results are shown in Figure 4. It is found
that the different simulation time scales show very different
patterns for the correlation function, the most obvious being an
uplifting of all of the correlation function curves, or rather, a
slower decay compared with the previous work.4,6,7,12 This point
can be approximately interpreted through the relation between
the positional fluctuation and the relaxation time τ of a given
atomic fluctuation correlation function, for which the initial
decay was approximated by an exponential function. Assuming
Figure 5. Atomic positional fluctuation autocorrelation functions for
10MeLeu-Cδ1 of CPA for the GLDBEM simulation. The results are
obtained by averaging on four different simulation periods in the 500ps simulation data: solid line, 100-140 ps; dotted line, 200-240 ps;
dashed line, 300-340 ps; and dash-dotted line, 400-440 ps.
that the atomic motions of a molecule obey a simple harmonic
Langevin equation, the relaxation time τ satisfies the realtion26,27
τ)γ
〈∆r2〉
3kBT
(30)
where T denotes the temperature, kB is Boltzmann’s constant,
〈∆r2〉 represents the positional fluctuation, and γ is the friction
coefficient. From our simulation, it is found that the longer the
simulation time, the larger the positional fluctuation. Therefore,
the results of Figure 4 are consistent with eq 30, in which the
400-ps simulation gave the largest positional fluctuation.
Another aspect concerns the stability and reliability of the
results. The results for the correlation functions obtained from
longer simulations in the present work are more stable than those
of the previous work mentioned above.
Here, it is necessary to provide a point of clarification that
the closeness between the two curves for the 200-ps and 400ps simulations in Figure 4 does not imply that there is an upper
limit beyond which longer simulations will provide essentially
the same positional fluctuation. According to our observations,
the curves for different atoms show significantly different
patterns; the curve for the 200-ps simulation is not always so
close to the curve for the 400-ps simulation. However, they do
have a common tendency: the longer the simulation time, the
more stable the correlation function curves.
Moreover, when averaged over the same simulation length
(40 ps) but at different time sections of the GLDBEM simulation, the correlation functions also show large differences (see
Figure 5). This means that the results for analyzing the
correlation function using the data from short simulations are
less reliable.
Conclusions
In this work, an efficient procedure for combining the BEM
with the GLD simulation technique has been described. The
two extra mean forces of the hydration force on a charge and
the surface pressure of the solvent determined by the BEM and
the memory effect of the friction kernel are together taken into
account in this method. The analysis results of our GLDBEM
simulation have been compared with those of MD, GLD, SD,
and SDBEM simulations. The results show that the extra mean
force and the friction memory effect of the solvent in the
simulations can increase the molecular flexibility and reduce
Hydration Force and Friction Memory Effect
the total number of intramolecular H-bonds. It was also found
that our GLDBEM simulation shows an obvious improvement
in dynamical properties compared with conventional SD and
SDBEM simulations. In addition, the results for the correlation
function obtained from 400-ps simulations in the present work
are much more stable than those from the previous 40-60-ps
simulations.
We would like to mention that some aspects of the present
method could be improved in future work: (1) The calculation
speed and precision of the algorithm for solving the PoissonBoltzmann equations using the BEM and of the numerical
integration scheme should be increased. In our work, the
simulation CPU time for each picosecond for the CPA molecule
on PII233 are as follows: 0.003 h for SD, 0.005 h for GLD,
0.085 h for GLDBEM, and 0.077 h for MD (with 754 water
molecules). The computation of solvation forces with the BEM
method takes significantly more time than the computation of
solute-solvent forces when solvent molecules are explicitly
represented. A new algorithm should be explored to accelerate
the BEM method. (2) Many MD simulations have shown that
the solvent dynamical behavior is bimodal rather than exponential. Therefore, the general form of the GLD simulation
should be studied under different models for the friction kernel.
(3) The relaxation times of the friction kernels are dependent
on atomic species. Therefore, different relaxation times could
be selected for the different atoms. (4) The anisotropy of the
stochastic force on an atom could be considered because of the
anisotropical exposure to solvent molecules.
Acknowledgment. We thank Professor W. F. van Gunsteren
for kindly providing the GROMOS package. C.X.W. thanks the
Abdus Salam International Centre for Theoretical Physics
(ICTP), Trieste, Italy, because this work was done (in part) in
the framework of Associate Membership Program of the ICTP.
This work was supported in part by the Chinese Natural Science
Foundation (No. 19774501, No. 29992590-2) and the Beijing
Natural Science Foundation (No. 5992002).
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6883
References and Notes
(1) Karplus, M.; Petsko, G. A. Nature 1990, 347, 631.
(2) van Gunsteren, W. F.; Weiner, P. K.; Wilkinson, A. J. Computer
Simulation of Biomolecular Systems: Theoretical and Experimental Application; Leidon ESCOM: The Netherlands, 1993; Vol. 2.
(3) Smith, J. C.; Karplus, M. J. Am. Chem. Soc. 1992, 114, 801.
(4) Wang, C. X.; Wan, S. Z.; Xiang, Z. X.; Shi, Y. Y. J. Phys. Chem.
1997, 101, 230.
(5) Fraternali, F.; van Gunsteren, W. F. J. Mol. Biol. 1996, 256, 939.
(6) Lautz, J.; Kessler, H.; van Gunsteren, W. F.; Weber, H. P.; Wenger,
R. M. Biopolymers 1990, 29, 1669.
(7) Shi, Y. Y.; Wang, L.; van Gunsteren, W. F. Mol. Simul. 1988, 1,
369.
(8) Zauhar, R. J. J. Comput. Chem. 1991, 12, 575.
(9) Sharp, K. J. Comput. Chem. 1991, 12, 454.
(10) Gilson, M. K.; Davis, M. E.; Luty, B. A.; McCammon, J. A. J.
Phys. Chem. 1993, 97, 3591.
(11) Gilson, M. K.; McCammon, J. A.; Madura, J. D. J. Comput. Chem.
1995, 16, 1081.
(12) Wan, S. Z.; Wang, C. X.; Shi, Y. Y. Mol. Phys. 1998, 93, 901.
(13) Xiang, Z. X.; Huang, F. H.; Shi, Y. Y. J. Phys. Chem. 1994, 98,
12782.
(14) Kubo, R. Rep. Prog. Phys. 1966, 29, 255.
(15) Vesely, F. J. Mol. Phys. 1984, 53, 505.
(16) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Simul. 1988, 1,
173.
(17) Papoulis, A. Probability, Random Variables and Stochastic Processes; 2nd ed.; McGraw-Hill: New York, 1985.
(18) Yoon, B. J.; Lenhoff, A. M. J. Comput. Chem. 1990, 11, 1080.
(19) Zauhar, R. J.; Morgan, R. S. J. Mol. Biol. 1985, 186, 815.
(20) Juffer, A. H.; Botta, E. F. F.; van Keulen, B. A. M.; Ploeg, A. V.
D.; Berendsen, H. J. C. J. Comput. Phys. 1991, 97, 144.
(21) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular
Simulation (GROMOS) Library Manual, BIOMOS: Groningen, The
Netherlands, 1987.
(22) Loosli, H. R.; Kessler, H.; Oschkinat, H.; Weber, H. P.; Petchez,
T. J.; Widmer, A. HelV. Chim. Acta 1985, 68, 661.
(23) Daggett, V.; Levitt, M. Annu. ReV. Biophys. Biomol. Struct. 1993,
22, 353.
(24) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola,
A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684.
(25) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys.
1977, 23, 327.
(26) Chandrasekchar, S. ReV. Mod. Phys. 1943, 15, 1.
(27) Swaminathan, S.; Ichiye, T.; van Gunsteren, W. F.; Karplus, M.
Biochemistry 1982, 21, 5230.