Molecular Dynamics Simulations of a Polyalanine Octapeptide under Ewald Boundary
Conditions: Influence of Artificial Periodicity on Peptide Conformation
Wolfgang Weber,† Philippe H. Hu1nenberger,*,‡ and J. Andrew McCammon†,⊥
Department of Chemistry and Biochemistry and Department of Pharmacology,
UniVersity of California at San Diego, 9500 Gilman DriVe, La Jolla, California 92093-0365 and Laboratorium
fu¨r Physikalische Chemie, ETH Zentrum, UniVersita¨tstrasse 6, CH-8092 Zu¨rich, Switzerland
ReceiVed: October 25, 1999; In Final Form: January 25, 2000
Ewald and related mesh methods are nowadays routinely used in explicit-solvent simulations of solvated
biomolecules, although they impose an artificial periodicity in systems which are inherently nonperiodic. In
the present study, we investigate the consequences of this approximation for the conformational equilibrium
of a polyalanine octapeptide (with charged termini) in water. We report three explicit-solvent molecular
dynamics simulations of this peptide in cubic unit cells of edges L ) 2, 3, and 4 nm, using the particleparticle-particle-mesh
(P3M) method for handling electrostatic interactions. The initial configuration of the
peptide is R-helical. In the largest unit cell (L ) 4 nm), the helix unfolds quickly toward configurations with
shorter end-to-end distances. By contrast, in the two smaller unit cells (L ) 2 and 3 nm), the R-helix remains
stable during 2 ns. Backbone fluctuations are somewhat larger in the medium (L ) 3 nm) compared to the
smallest unit cell. These differences are rationalized using a continuum electrostatics analysis of configurations
from the simulations. These calculations show that the R-helical conformation is stabilized by artificial
periodicity relative to any other configuration sampled during the trajectories. This artificial stabilization is
larger for smaller unit cells, and is responsible for the absence of unfolding in the two smaller unit cells and
the reduced backbone fluctuations in the smallest unit cell. These results suggest that artificial periodicity
imposed by the use of infinite periodic (Ewald) boundary conditions in explicit-solvent simulations of
biomolecules may significantly perturb the potentials of mean force for conformational equilibria, and even
in some cases invert the relative stabilities of the folded and unfolded states.
Introduction
The treatment of long-range electrostatic interactions in
molecular simulations of liquids and solutions has been an area
of active research for a long time.1-6 The main reasons are that
(i) this treatment generally represents the computationally most
expensive part of a simulation and (ii) simulated observables
are often extremely sensitive to approximations made in this
treatment. During the past few years, the Ewald7 and related
mesh methods (P3M8 and PME10) to treat long-range electrostatic
interactions in molecular simulations have been implemented
in many of the widely used simulation packages (e.g.,
CHARMm, AMBER, NWChem). These methods are nowadays
routinely used in explicit-solvent simulations of biomolecules11-19
and may soon completely supersede the use of cutoff-based
methods, which rely on the truncation of long-range electrostatic
interactions. While the artifacts caused by cutoff truncation
have been extensively studied11-13,20-34 and shown to be severe
in many cases, possible artifacts linked with the use of the
Ewald and related mesh methods are just beginning to be
understood.6,35-42 Since the reliability of computer simulations
may be significantly impaired by errors committed in the
treatment of electrostatic interactions, it is of importance to
assess also these methods with respect to the nature and
magnitude of artifacts they may induce in various simulated
observables.
Ewald and related mesh methods permit, in principle, the
exact evaluation of electrostatic interaction energies and forces
in systems with infinite periodic (Ewald) boundary conditions.
However, whereas infinite translational symmetry is probably
a good approximation for the simulation of crystals, it may prove
a bad approximation for the simulation of liquids and solutions.6,41,42
In particular, it can be expected that imposing
artificial periodicity on liquids and solutions will represent a
negligible perturbation when the unit cell is very large, but will
drastically affect the properties of the system when small unit
cells are used. The minimal size of the unit cell necessary to
ensure a negligible periodicity-induced perturbation will depend
on the nature of the system and on the observables of interest.
A number of studies suggested that artifacts linked with Ewald
boundary conditions for the simulation of solutions are negligible.
Among these, studies of the dependence of the potential
energy on the orientation of small biomolecules within the unit
cell indicated that Ewald rotational artifacts are small when a
solvent of high permittivity is considered.43,44 However, this
conclusion may be restricted to solutes with a small solventexcluded
volume compared to the size of the unit cell.42,45 A
comparison of three 1.5 ns explicit-solvent simulations of a DNA
dodecamer using the PME method and different numbers of
water molecules46 did not evidence a significant unit-cell size
dependence of the monitored properties (root-mean-square
deviation from the crystallographic structure and DNA curvature).
This result may indicate the absence of artifacts linked
with Ewald boundary conditions in this specific system, but this
absence of artifacts may also simply be caused by a poor
convergence or weak sensitivity of the selected observables.
* Corresponding author. E-mail: phil@igc.phys.chem.ethz.ch.
† Department of Chemistry and Biochemistry, University of California
at San Diego.
‡ ETH Zentrum.
⊥ Department of Pharmacology, University of California at San Diego.
3668 J. Phys. Chem. B 2000, 104, 3668-3675
10.1021/jp9937757 CCC: $19.00 © 2000 American Chemical Society
Published on Web 03/22/2000
Recently, we proposed a new method based on continuum
electrostatics for assessing the nature and magnitude of periodicity-induced
artifacts in explicit-solvent simulations under Ewald
boundary conditions.6,41,42 The principle of these calculations
is the following. Ideally, one would like to perform a direct
comparison between explicit-solvent simulations under nonperiodic
and Ewald boundary conditions. This is impossible
since an explicit-solvent system under nonperiodic boundary
condition would be prohibitively large. However, a comparison
between nonperiodic and Ewald boundary conditions becomes
possible when an implicit representation of the solvent is used
instead, e.g., by employing continuum electrostatics. Assuming
that observations made using continuum electrostatics are (at
least qualitatively) relevant for explicit-solvent simulations, this
approach allows to investigate specifically the perturbation
induced by artificial periodicity in simulations under Ewald
boundary conditions.
In contrast to the previously mentioned studies,43,44,46 the
continuum electrostatics analysis of potentials of mean force
for model conformational transitions of biomolecules indicated
that artificial periodicity may significantly affect the result of
biomolecular simulations.6,42 The continuum electrostatics analysis
of configurations from a nanosecond molecular dynamics
(MD) simulation (using the P3M method) of the protein Sac7d
at 550 K47 suggested that the anomalous stability of the protein
during this simulation is (at least in part) a consequence of
artificial periodicity.42
In the present work, we investigate whether the qualitative
conclusions reached using the continuum electrostatics study
of a model conformational transition can be confirmed by
explicit-solvent MD simulations. Here, we focus on the case of
a polyalanine octapeptide with charged termini, for which strong
periodicity-induced artifacts were observed.42 We report three
explicit-solvent MD simulations of this peptide using the P3M
method for handling electrostatic interactions and employing
three different unit cell sizes (cubes of edges L ) 2, 3, and 4
nm). The effect of artificial periodicity on the conformational
behavior of the peptide is analyzed in detail. It is then shown
how the outcome of the simulations can be rationalized with
the help of the continuum electrostatics analysis.
Theory
Using a continuum representation of the solvent, the free
energy of a solute-solvent system may be partitioned into a
sum of a nonpolar contribution, corresponding to the work
required for creating the solute cavity in the solvent in the
absence of solute charges, and an electrostatic contribution,
corresponding to the work required for reversibly charging the
solute atoms. Since the nonpolar contribution is dominated by
short-range effects, we assume that the periodicity-induced
perturbation of this contribution is negligible. Therefore, we
focus on the perturbation imposed by artificial periodicity on
the electrostatic contribution alone. This perturbation can be
quantified as
where [∆Gel]NP and [∆Gel]P are the electrostatic free energies
of the solute-solvent system under nonperiodic and Ewald
boundary conditions, respectively. As described elsewhere,6,41,42
these quantities can be partitioned into two and three terms,
respectively, as
and
In practice, the solute is described as a low-dielectric cavity of
relative permittivity i immersed in a solvent of higher permittivity
s. The Coulomb terms ECb represent interactions between
the partial charges on the solute atoms in a homogeneous
medium of permittivity i. The solvation terms ∆Gsolv correspond
to the work required for changing the external permittivity from
i to s. Finally, the self-energy term ∆Gself (Ewald boundary
conditions) corresponds to the sum of the self-energies of the
charges plus homogeneous neutralizing background density
(C + B charges41). A more detailed discussion of these different
terms, together with numerical methods for their evaluation, is
given elsewhere.41
Computational Details
Here, we investigate the unfolding equilibrium of a polyalanine
octapeptide with charged N- and C-termini.
In a first step, we use the continuum electrostatics method to
analyze the periodicity-induced perturbation of the electrostatic
free energy (eq 1) along a model unfolding pathway. Peptide
conformations along this pathway are generated using the
dimensionless parameter λ, which defines the backbone φ and
ψ angles of all residues through
Other internal coordinates (bond lengths, bond angles, and
improper dihedral angles) are fixed to standard values from the
GROMOS96 force field.48 The values λ ) 0 and λ ) 1
correspond to the canonical R-helix conformation and to the
extended (all-trans) conformation, respectively. For all conformations,
the three principal axes of the peptide are aligned with
an edge of the unit cell. The periodicity-induced perturbation
is calculated for cubic unit cells of edges L ) 2, 3, and 4 nm.
The parameters for these calculations are the same as for the
analysis of the configurations from the MD simulations (see
below).
In a second step, we report the results of explicit-solvent MD
simulations of the peptide. These simulations were performed
using the GROMOS96 force field49 and a modified version
of the GROMOS96 simulation program48,50 incorporating the
P3M method8 for handling electrostatic interactions. The initial
configuration was taken to be the canonical R-helix (eq 4 with
λ ) 0). The solute was immersed in three different cubic unit
cells containing 234, 866, and 2099 SPC water molecules51 so
that the edges of the corresponding unit cells were L ) 2, 3,
and 4 nm, respectively, with a solvent density of 1.0 g‚cm-3.
In the starting configuration, the minimal distances between any
solute atom and the nearest cell wall were 0.4, 1.4, and 2.4 nm,
respectively. The three simulations were performed under Ewald
boundary conditions using the P3M method,8 with a 1.0 nm
cutoff distance for the real space contribution, a grid spacing
of 0.05 nm, and a Gaussian charge-shaping function of width
0.3 nm. A 1.0 nm cutoff was used for the Lennard-Jones
interactions. The nonbonded pair list was updated every 10 steps.
The simulations were carried out in the NVT ensemble using a
Berendsen thermostat52 with separate solute and solvent coupling,
a reference temperature of 300 K, and a relaxation time
of 0.1 ps. Bond lengths were constrained using SHAKE53 with
a relative tolerance of 10-4. A 2 fs time step was used to
integrate the equations of motion, and coordinates were saved
∆∆Gel ) [∆Gel]P - [∆Gel]NP (1)
[∆Gel]NP ) [ECb]NP + [∆Gsolv]NP (2)
[∆Gel]P ) [ECb]P + [∆Gself]P + [∆Gsolv]P (3)
φ(λ) ) -57°(1 - λ) - 180° λ and
ψ(λ) ) -47°(1 - λ) - 180° λ (4)
Molecular Dynamics Simulations of an Octapeptide J. Phys. Chem. B, Vol. 104, No. 15, 2000 3669
every 0.1 ps for analysis. For equilibration, the system was
heated to 300 K in steps of 100 K, each of 5 ps duration.
Positional restraints of decreasing magnitude were applied
simultaneously with the heating process. After equilibration, data
were collected for 2 ns (L ) 2 and 3 nm) or 1 ns (L ) 4 nm).
For the analysis of the trajectories, we monitored the CR rootmean-square
atomic positional deviations (drms) with respect to
the structure of the canonical R-helix, the CR root-mean-square
atomic positional fluctuations (frms), the hydrogen-bonding
pattern, and the end-to-end distance (measured between the
N-terminal ammonium nitrogen and the C-terminal carboxlyate
carbon). For the analysis of the hydrogen-bonding pattern, a
hydrogen bond was assumed to be formed if the distance
between the hydrogen and acceptor atoms was smaller than 0.25
nm and the angle between donor, hydrogen, and acceptor was
larger than 135°. Hydrogen bonds corresponding to the R-helical
configuration were numbered 1 through 4 according to the
following definitions (residue numbers are indicated in parentheses):
1 ) NH(5)fCO(1), 2 ) NH(6)fCO(2), 3 ) NH(7)fCO(3),
and 4 ) NH(8)fCO(4).
Finally, configurations sampled at 10 ps intervals during the
three simulations were analyzed using the continuum electrostatics
method to extract the periodicity-induced perturbation
of the electrostatic free energy ∆∆Gel (eq 1). To this purpose,
the electrostatic free energies under nonperiodic (eq 2) and
Ewald (eq 3) boundary conditions were calculated as described
previously41,42 using a modified version of the UHBD program54
with a grid spacing of 0.05 nm. The calculations under Ewald
boundary conditions were performed using cubic unit cells of
edges L ) 2, 3, and 4 nm, as appropriate. For the calculations
under nonperiodic boundary conditions, the edge lengths were
augmented by a distance of 0.2 nm with respect to the value
used in the simulation, in order to increase the distances between
the solute atoms and the boundary of the grid.42 Solute cavities
were defined as the contact and reentrant surface obtained by
rolling a probe of 0.14 nm radius over the peptide. For
consistency, the atomic radii and charges were taken from the
GROMOS96 force field.48,49 The solute and solvent permittivities
were set to i ) 1 and s ) 78, respectively.
Results and Discussion
First, we briefly discuss the results obtained from the analysis
of the model unfolding pathway defined by eq 4. The periodicity-induced
perturbation of the electrostatic free energy ∆∆Gel
(eq 1) along this pathway is displayed in Figure 1 as a function
of the end-to-end distance. The quantity ∆∆Gel is always
negative and its absolute value increases with increasing endto-end
distance. Thus, artificial periodicity tends to destabilize
the R-helical conformation relative to more extended structures.
As shown by the different slopes of the curves corresponding
to increasing edge lengths of the cubic unit cell (L ) 2 nm to
L ) 4 nm), this effect becomes less important when the unit
cell size is increased. This is not surprising since finite-size
effects should vanish completely in the limit of an inifitely large
unit cell.
As illustrated in Figure 2 for L ) 4 nm, the trends in ∆∆Gel
result from the balance between two opposing effects.42 On the
one hand, the Coulomb contributions [ECb]NP and [ECb]P +
[∆Gself]P (Figure 2a) decrease with increasing end-to-end
distance, because the unfavorable interaction between the
backbone dipole and the charged termini decreases in magnitude
when the peptide unfolds. The shift in the Coulomb contribution
caused by artificial periodicity is negative and its absolute value
increases in magnitude with increasing end-to-end distance. This
shift is due to the favorable interaction between the charged
termini of the oligopeptide and those of its periodic copies in
adjacent unit cells. On the other hand, the solvation contributions
[∆Gsolv]NP and [∆Gsolv]P (Figure 2b) are negative and their
absolute values decrease with increasing end-to-end distance.
The shift in the solvation contribution caused by artificial
periodicity is positive and increases in magnitude with increasing
Figure 1. Periodicity-induced perturbation of the electrostatic free
energy ∆∆Gel (eq 1), for a polyalanine octapeptide (charged termini),
along the model unfolding pathway defined by eq 4. The end-to-end
distance is measured between the N-terminal ammonium nitrogen and
the C-terminal carboxlyate carbon. Curves are displayed for cubic unit
cells of edges L ) 2, 3, and 4 nm. The calculations were performed as
described in ref 42, using the GROMOS96 charges and radii.48,49 These
curves differ from those reported in ref 42 and result from more accurate
calculations.
Figure 2. Coulomb (a) and solvation (b) contributions to the
electrostatic free energy (c) calculated under nonperiodic (NP) or Ewald
(P) boundary conditions (see eqs 2 and 3) for a polyalanine octapeptide
(charged termini), along the model unfolding pathway defined by eq
4. Curves are displayed for a cubic unit cell of edge L ) 4 nm.
3670 J. Phys. Chem. B, Vol. 104, No. 15, 2000 Weber et al.
end-to-end distance. This shift may be viewed as a consequence
of the perturbation of the solvent by the periodic copies of the
oligopeptide, which renders the solvent less available for the
solvation of the oligopeptide in the central unit cell. This effect
increases in magnitude with the size of the overall dipole
moment of the peptide, i.e., with the end-to-end distance. The
sum of the Coulomb and solvation contributions, i.e., the total
electrostatic contribution to the potential of mean force, is shown
in Figure 2c for nonperiodic and Ewald boundary conditions.
This total electrostatic contribution strongly favors extended
conformations relative to the R-helix. Note, however, that the
overall potential of mean force depends on additional covalent,
nonpolar, and entropic contributions, so that no direct conclusion
can be drawn on which conformations will dominate at
equilibrium. The large periodicity-induced shifts in the Coulomb
and solvation contributions (Figure 2a,b) cancel each other to
a large extent, resulting in a much smaller periodicity-induced
shift of the overall electrostatic free energy (Figure 2c). The
periodicity-induced shift in the Coulomb contribution (favoring
extended structures) still slightly dominates over the corresponding
shift in the solvation contribution (favoring the
R-helix), resulting in the trends observed in Figure 1.
To check whether these qualitative observations reached using
continuum electrostatics applied to a model unfolding pathway
can be verified in explicit-solvent MD simulations, we now
discuss the results of three explicit-solvent MD simulations of
the octapeptide in cubic unit cells of edges L ) 2, 3, and 4 nm,
and using the P3M method for handling electrostatic interactions.
The time evolutions of the CR root-mean-square atomic positional
deviation (drms, lower panels) and of the hydrogen-bonding
pattern (upper panels) are displayed in Figure 3 for the three
simulations. For the L ) 2 nm simulation, drms remains relatively
low (about 0.1 nm) and fluctuates little over the 2 ns simulation
time. The same observation applies to the L ) 3 nm simulation,
except for two excursions to higher values (about 0.25 nm) in
the intervals 0.35-0.45 and 1.50-1.60 ns. By contrast, the value
of drms during the L ) 4 nm simulation increases rapidly to
about 0.4 nm within the first 0.6 ns of the simulation. A
qualitative difference can also be observed between the hydrogenbonding
patterns corresponding to the three simulations. While
the C-terminal hydrogen bond 4 is disrupted very early in all
cases, the three other R-helical hydrogen bonds are only
preserved during the L ) 2 and 3 nm simulations. For L ) 2
nm, hydrogen bonds 1, 2, and 3 are present for 90%, 88%, and
82% of the configurations, respectively. Hydrogen bonds
characteristic of a 310-helix (not shown in the figure), such as
NH(5)fCO(2), NH(6)fCO(3), and NH(7)fCO(4), also occur
during this simulation, but are formed in less than 5% of the
configurations. For L ) 3 nm, hydrogen bonds 1, 2, and 3 are
present for 75%, 76%, and 67% of the configurations, respectively.
Transient disruptions of these three hydrogen bonds
coincide with the two peaks observed in drms around 0.40 and
1.55 ns. Here, 310-helix hydrogen bonds are slightly more
frequent (5-10% of the configurations). In addition, the
hydrogen bonds NH(6)fCO(1) and NH(7)fCO(2) replace the
R-helical ones around 0.40 and 1.55 ns. For L ) 4 nm, the
initial R-helical hydrogen bonds are disrupted after about 0.06
ns (bonds 2 and 3) and 0.73 ns (bond 1). No other type of
hydrogen bond occurs for more than 5% of the configurations,
except for the hydrogen bond NH(7)fCO(1) present from 0.75
ns until the end of the simulation.
The end-to-end distance of the peptide is displayed as a
function of time in Figure 4 (upper panels). For the simulations
in the two smallest unit cells, the end-to-end distance remains
close to the initial value of 1.27 nm (canonical R-helix), as could
be expected from Figure 3. By contrast, for L ) 4 nm, an initial
increase of the end-to-end distance during the first 0.1 ns is
followed by a rapid decrease to lower values in the range 0.3-
1.1 nm, compatible with the existence of the NH(7)fCO(1)
hydrogen bond. The periodicity-induced perturbation of the
electrostatic free energy ∆∆Gel (eq 1) is also displayed in Figure
4 (lower panels). Three observations can be made here: (i)
∆∆Gel is always negative; (ii) it decreases in magnitude with
increasing unit cell size; (iii) it is strongly anti-correlated with
the end-to-end distance. The latter point is particularly evident
in the L ) 4 nm simulation where the octapeptide undergoes a
significant conformational change. These three observations are
in perfect agreement with the results obtained from the study
of the model unfolding pathway (Figure 1).
Ten superimposed configurations sampled at 100 ps intervals
during each of the three simulations are shown in Figure 5. The
corresponding CR root-mean-square atomic positional fluctuations
(frms) are displayed in Figure 6. For L ) 2 nm, as could
be anticipated from Figures 3a and 4a, the peptide remains very
stable in the R-helical conformation. The corresponding backbone
fluctuations are small (0.04-0.06 nm for the CR of residues
1-7). For L ) 3 nm, the peptide also remains stable in the
R-helical conformation. However, the backbone fluctuations are
somewhat higher (0.07-0.09 nm for the CR of residues 1-7)
and transient disruption of R-helical hydrogen bonds are
observed (Figure 3b). By contrast, for L ) 4 nm, the peptide is
completely unfolded after about 0.7 ns (Figure 3c), and the
backbone fluctuations are significantly higher than in the two
previous cases. Since the configurations from the last 0.3 ns of
this simulation lack any hydrogen bond except for NH(7)fCOFigure
3. CR root-mean-square atomic positional deviation (drms) from
the canonical R-helix (lower panels) and hydrogen-bonding pattern
(upper panels), from the MD simulations with cubic unit cells of edges
L ) 2 (a), 3 (b), and 4 nm (c). The hydrogen bond numbering scheme
is (residue numbers are given in parentheses): 1 ) NH(5)fCO(1), 2
) NH(6)fCO(2), 3 ) NH(7)fCO(3), and 4 ) NH(8)fCO(4).
Molecular Dynamics Simulations of an Octapeptide J. Phys. Chem. B, Vol. 104, No. 15, 2000 3671
(1) and have a short end-to-end distance (Figure 4c) so that the
two charged termini are almost in direct contact (Figure 5c),
we will refer to these as “circle” conformations. Note that in
contrast to the simplistic unfolding pathway used in the initial
continuum electrostatics study (Figure 1), where unfolding was
interpreted as an increase in the end-to-end distance, unfolding
in the L ) 4 nm simulation leads to circle conformations with
shorter end-to-end distances compared with the R-helix. Except
for this important difference, the results of the explicit-solvent
MD simulations and of the continuum calculations on the model
unfolding pathway are in very good qualitative agreement, as
discussed below.
The cause for the different conformational behavior of the
polypeptide chain in the three simulations can be easily
understood by considering Figure 7. This figure shows, for the
conformations sampled at 10 ps intervals during the three
simulations, the correlation between the periodicity-induced
perturbation of the electrostatic free energy ∆∆Gel (calculated
using the continuum method) and the end-to-end distance. The
points representative for the canonical R-helix and for a typical
circle configuration (sampled at 0.8 ns during the L ) 4 nm
simulation) are also represented. Note the similarity between
Figures 1 and 7 in terms of the slopes of the correlations for
the three different unit cell sizes. As discussed above (see Figure
4), the absolute value of ∆∆Gel decreases with both a decreasing
end-to-end distance and an increasing unit cell size. The
periodicity-induced destabilization of the circle configuration
relative to the canonical R-helix evaluates to 2.8 (L ) 2 nm),
0.7 (L ) 3 nm), and 0.3 kJ/mol (L ) 4 nm). Configurations
which are penalized by a perturbation of the order of (or greater
than) kBT ) 2.5 kJ/mol at 300 K are essentially inaccessible
during the MD simulation. For L ) 2 nm, the free energy
corresponding to configurations with shorter end-to-end distances
(compared to the R-helix) is strongly affected by artificial
periodicity. Since the circle configurations are penalized by
about 2.8 kJ/mol (greater than kBT ), they are essentially
inaccessible during the MD simulation. For L ) 3 nm, ∆∆Gel
as well as the slope of the correlation line are smaller in
magnitude, and the circle configurations are penalized by only
about 0.7 kJ/mol. Consequently, configurations with shorter endto-end
distances become more accessible, although not to the
extent of permitting unfolding within the 2 ns simulation. On
the other hand, fluctuations in the peptide backbone become
more important (Figure 6). For L ) 4 nm, the periodicityinduced
perturbation becomes small enough to permit unfolding
toward the (thermodynamically more stable) circle configuration.
The behavior of the peptide in this latter simulation is probably
the most relevant for comparison with experiment, since the
perturbation due to artificial periodicity is the smallest in this
case. Indeed, it is known experimentally that oligopeptides of
such short lengths do not adopt a stable R-helical conformation.55
Note finally that the quality of the correlation between ∆∆Gel
and the end-to-end distance improves in the sequence L ) 2
nm (correlation coefficient r ) -0.52 ), L ) 3 nm (r ) -0.73),
and L ) 4 nm (r ) -0.94). This can be understood as follows.
In the largest box, periodicity-induced artifacts are caused
primarily by the interaction between the peptide in the central
unit cell and its periodic copies in the two adjacent cells along
the peptide axis, and thus are well correlated with the end-toend
distance. On the other hand, in the smallest box, the
interaction between the peptide and its periodic copies in the
other four adjacent cells becomes more important (because these
copies are now only about 2 nm away), and the correlation with
the end-to-end distance becomes worse.
The above explanation for the different conformational
behaviors of the peptide during the three simulations rests on
the hypothesis that the circle configuration is thermodynamically
favored over the R-helical one. However, it is difficult to
distinguish whether periodicity-induced artifacts cause these
different behaviors for thermodynamic or kinetic reasons. If the
thermodynamic explanation holds, the free energy difference
between circle and helical conformations must be small enough
that the periodicity-induced perturbation in the smallest unit cells
is suffient to inverse the relative stabilities of the two conformations.
If the kinetic explanation holds, the periodicity-induced
perturbation in the smallest unit cells only slows down sampling
in the direction of the favored circle configuration, so that this
configuration is not observed during a nanosecond simulation.
Considering the magnitude of the periodicity-induced perturbation
(Figure 7), the thermodynamic explanation seems more
likely in the L ) 2 nm case, whereas the kinetic explanation is
probably correct for L ) 3 nm. This point could be addressed
in more detail by performing longer simulations. Note finally
that the peptide is never observed to unfold toward extended
configurations, even though these would be significantly
stabilized by artificial periodicity in the smallest boxes. This
probably simply reflects the intrinsic thermodynamic instability
of these configurations in the case of an octapeptide with
charged termini.
Another possible explanation for the absence of unfolding
in the two smaller unit cells (and especially for L ) 2 nm) could
be the lack of space, i.e., the presence of short-range packing
contacts between the peptide in the central unit cell and its
Figure 4. Periodicity-induced perturbation of the electrostatic free
energy ∆∆Gel (lower panels), according to eq 1, and end-to-end distance
e (upper panels), from MD simulations with cubic unit cells of edges
L ) 2 (a), 3 (b), and 4 nm (c). Dashed lines represent the ∆∆Gel and
e values characteristic of the R-helix: e ) 1.27 nm and ∆∆Gel ) -3.46
(L ) 2 nm), -0.85 (L ) 3 nm), and -0.33 kJ/mol (L ) 4 nm). For the
calculation of ∆∆Gel using continuum electrostatics, the configurations
were translated but not rotated, and the same unit cell sizes were used
as during the corresponding simulation.
3672 J. Phys. Chem. B, Vol. 104, No. 15, 2000 Weber et al.
periodic copies. However, this explanation can easily be ruled
out. The maximal values of the end-to-end distance in the three
simulations are (see Figure 4) 1.69 (L ) 2 nm), 1.78 (L ) 3
nm), and 1.85 nm (L ) 4 nm), which is always at least 0.3 nm
smaller than the edge length of the unit cell, and at least 1.7
nm smaller than the diagonal of the cell ( 3L). In fact, for the
L ) 2 nm simulation, the distance between the N-terminus of
the peptide in the central unit cell and the C-terminus of any of
its periodic copies is always larger than 0.5 nm, and larger than
0.75 nm in 90% of the configurations. Furthermore, the fact
that the unfolding observed in the L ) 4 nm simulation proceeds
with a maximal end-to-end distance of 1.85 nm indicates that
this pathway would be feasible in the L ) 2 nm unit cell without
involving any packing contacts. Thus, the reason for the absence
of unfolding in the two smaller unit cells is not the presence of
packing contacts within the periodic system, but rather the effect
of the periodicity-induced perturbation of the electrostatic
interactions.
Conclusion
In the present study we investigated artifacts linked with the
use of Ewald boundary conditions in explicit-solvent simulations
of a polyalanine octapeptide with charged termini. This was
done by comparing the results of three explicit-solvent MD
simulations of the peptide using the P3M method and cubic
unit cells of different sizes (L ) 2, 3, and 4 nm). These
Figure 5. Ten configurations (superimposed on the CR atoms) taken at 100 ps intervals from the three MD simulations with cubic unit cells of
edges L ) 2 (a), 3 (b), and 4 nm (c).
Figure 6. CR root-mean-square atomic positional fluctuations (frms)
corresponding to the eight residues of the octapeptide, calculated from
the first 1 ns of the three MD simulations with cubic unit cells of edges
L ) 2, 3, and 4 nm.
Figure 7. Periodicity-induced perturbation of the electrostatic free
energy ∆∆Gel (eq 1) calculated from configurations sampled at 10 ps
interval during the three MD simulations with cubic unit cells of edges
L ) 2, 3, and 4 nm. Values corresponding to the canonical R-helix
(×) and to a model circle configuration (+) are also represented (from
bottom to top L ) 2, 3, and 4 nm). The model circle configuration is
the configuration sampled at 0.8 ns during the MD simulation with L
) 4 nm. See also Figure 4.
Molecular Dynamics Simulations of an Octapeptide J. Phys. Chem. B, Vol. 104, No. 15, 2000 3673
simulations, initiated using the R-helical conformation, were
analyzed in terms of the conformational behavior of the peptide
chain.
Experimentally, peptides of such a small size are known not
to adopt stable helical conformations in water at room temperature.55
Our results show that indeed, under the most realistic
simulation conditions (i.e., in the largest unit cell of edge L )
4 nm), the R-helix unfolds quickly. The peptide was found to
adopt circle configurations in which the two charged termini
are in close proximity and no hydrogen bonds are present except
for NH(7)fCO(1). By contrast, during the simulations in the
two smaller unit cells (L ) 2 and 3 nm), which are more
significantly affected by finite-size effects, the R-helix remained
stable over the simulation times of 2 ns. The two latter
simulations differed in terms of the flexibility of the peptide
backbone, which was found to be more pronounced in the L )
3 nm simulation.
The different conformational behavior of the peptide observed
in these three MD simulations was rationalized using a
continuum electrostatic analysis of configurations sampled
during the simulations. Such calculations6,41,42 permit estimation
of the periodicity-induced perturbation of the electrostatic free
energy ∆∆Gel (eq 1) by comparing the electrostatic free energies
of the solute-solvent system under nonperiodic and Ewald
boundary conditions. The analysis of ∆∆Gel for the three
different simulations showed that the canonical R-helix is
stabilized by artificial periodicity relative to all other configurations
sampled during the trajectories. In particular, the R-helix
is stabilized relative to the circle configurations encountered in
the L ) 4 nm simulation after unfolding. This periodicityinduced
destabilization was estimated to be about 2.8 (L ) 2
nm), 0.7 (L ) 3 nm), and 0.3 kJ/mol (L ) 4 nm). This
perturbation appears sufficient to prevent (either thermodynamically
or kinetically) the unfolding in the L ) 2 and 3 nm
simulations. Furthermore, it significantly reduces the backbone
fluctuations in the L ) 2 nm simulation. Only for L ) 4 nm is
the perturbation sufficiently small so that unfolding readily
occurs.
These results suggest that artificial periodicity imposed by
the use of Ewald boundary conditions in explicit-solvent
simulations may significantly perturb the potentials of mean
force for conformational equilibria of solvated biomolecules.
In the case of the polyalanine octapeptide and in the two smaller
unit cells (L ) 2 and 3 nm), the perturbation may even be strong
enough to inverse the relative stabilities of the circle and
R-helical conformations. Considering the L ) 4 nm system as
a reference, it appears that periodicity-induced artifacts can be
made essentially negligible by sufficiently solvating the peptide
during the setup of the simulation. For this specific system, we
would recommend the inclusion of at least three solvation layers
(about 0.85 nm) between the peptide termini and the nearest
unit-cell walls. In fact, such a requirement is not unrealistically
drastic for the simulation of typical biomolecules.
Considering the results of the present study, we suggest that
explicit-solvent simulations of biomolecular systems under
Ewald boundary conditions should be systematically examined
for periodicity-induced artifacts. The continuum electrostatics
analysis proposed previously6,41,42 appears to be an appropriate
(qualitative) tool for this purpose. From this and previous
studies,41,42 periodicity artifacts are likely to be important (kBT
or larger) for systems involving (i) a solvent of low dielectric
permittivity, (ii) few solvation shells around the solute, and (iii)
a solute bearing a net charge or a large dipole. In particular,
the periodicity artifacts observed in the present study would
probably be negligible in the case of an octapeptide with neutral
termini.42
However, in the long run, we also hope that the realization
that periodicity-induced artifacts are sometimes far from
negligible will trigger the development of more accurate latticesum
methods for handling long-range interactions in simulations.
Acknowledgment. P.H.H. acknowledges support from the
Swiss National Science Foundation and the Human Frontier
Science Program. Additional support was provided by NIH and
NSF. This work was also partially supported by the NSF
Supercomputer Centers NRAC program.
References and Notes
(1) Harvey, S. C. Proteins: Struct. Funct. Genet. 1989, 5, 78.
(2) Davis, M.E.; McCammon, J.A. Chem. ReV. 1990, 90, 509.
(3) Smith, P.E. In Computer simulation of biomolecular systems,
theoretical and experimental applications; van Gunsteren, W. F., Weiner,
P. K., Wilkinson, A. J., Eds.; ESCOM: Leiden, The Netherlands, 1993;
Vol. II, pp 182-212.
(4) Levy, R. M.; Gallicchio, E. Annu. ReV. Phys. Chem. 1998, 49, 531.
(5) Sagui, C.; Darden, T. A. Annu. ReV. Biophys. Biomol. Struct. 1999,
28, 155.
(6) Hu¨nenberger, P. H. In Simulation and theory of electrostatic
interactions in solution: Computational chemistry, biophysics, and aqueous
solutions; Pratt, L. R., Hummer, G., Eds.; American Institute of Physics:
New York, 1999; pp 17-83.
(7) Ewald, P. P. Ann. Phys. 1921, 64, 253.
(8) Hockney, R. W.; Eastwood, J. W. Computer simulation using
particles; Institute of Physics Publishing: Bristol, 1988.
(9) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089.
(10) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.;
Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577.
(11) Smith, P. E.; Pettitt, B. M. J. Chem. Phys. 1991, 95, 8430.
(12) Schreiber, H.; Steinhauser, O. Biochemistry 1992, 31, 5856.
(13) Schreiber, H.; Steinhauser, O. Chem. Phys. 1992, 168, 75.
(14) York, D. M.; Darden, T. A.; Pedersen, L. G. J. Chem. Phys. 1993,
99, 8345.
(15) Fox, T.; Kollman, P. A. Proteins: Struct. Funct. Genet. 1996, 25,
315.
(16) Mohan, V.; Smith, P. E.; Pettitt, B. M. J. Phys. Chem. 1993, 97,
12984.
(17) Cheatham III, T. E.; Kollman, P. A. J. Mol. Biol. 1996, 259, 434.
(18) Young, M. A.; Ravishanker, G.; Beveridge, D. L. Biophys. J. 1997,
73, 2313.
(19) Auffinger, P.; Westhof, E. Curr. Opin. Struct. Biol. 1998, 8, 227.
(20) Neumann, M. Mol. Phys. 1983, 50, 841.
(21) Neumann, M.; Steinhauser, O.; Pawley, G. S. Mol. Phys. 1984,
52, 97.
(22) Wood, R. H. J. Chem. Phys. 1995, 103, 6177.
(23) Brooks III, C. L.; Pettitt, B. M.; Karplus, M. J. Chem. Phys. 1985,
83, 5897.
(24) Baker, N. A.; Hu¨nenberger, P. H.; McCammon, J. A. J. Chem.
Phys. 1999, 110, 10679.
(25) Brooks III, C. L. J. Chem. Phys. 1987, 86, 5156.
(26) Madura, J. D.; Pettitt, B. M. Chem. Phys. Lett. 1988, 150, 105.
(27) Straatsma, T. P.; Berendsen, H. J. C. J. Chem. Phys. 1989, 89,
5876.
(28) Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1986, 84, 5836.
(29) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1990, 94, 4303.
(30) Gua`rdia, E.; Rey, R.; Padro´, J. A. J. Chem. Phys. 1991, 95, 2823.
(31) Dang, L. X.; Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1992, 96,
4046.
(32) Hummer, G.; Soumpasis, D. M.; Neumann, M. Mol. Phys. 1993,
81, 1155.
(33) Bader, J. S.; Chandler, D. J. Phys. Chem. 1992, 96, 6423.
(34) Schreiber, H.; Steinhauser, O. J. Mol. Biol. 1992, 228, 909.
(35) Figueirido, F.; Del Buono, G. S.; Levy, R. M. J. Chem. Phys. 1995,
103, 6133.
(36) Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Phys. Chem. 1996, 100,
1206.
(37) Figueirido, F.; Del Buono, G. S.; Levy, R. M. J. Phys. Chem. B
1997, 101, 5622.
(38) Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Chem. Phys. 1997, 107,
9275.
3674 J. Phys. Chem. B, Vol. 104, No. 15, 2000 Weber et al.
(39) Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Phys. Chem. A 1998,
102, 7885.
(40) Sakane, S.; Ashbaugh, H. S.; Wood, R. H. J. Phys. Chem. B 1998,
102, 5673.
(41) Hu¨nenberger, P. H.; McCammon, J. A. J. Chem. Phys. 1999, 110,
1856.
(42) Hu¨nenberger, P. H.; McCammon, J. A. Biophys. Chem. 1999, 78,
69.
(43) Smith, P. E.; Pettitt, B. M. J. Chem. Phys. 1996, 105, 4289.
(44) Smith, P. E.; Blatt, H. D.; Pettitt, B. M. J. Phys. Chem. B 1997,
101, 3886.
(45) Luty, B. A.; van Gunsteren, W. F. J. Phys. Chem. 1996, 100, 2581.
(46) Norberto de Souza, O.; Ornstein, R. L. Biophys. J. 1997, 72, 2395.
(47) de Bakker, P. I. W.; Hu¨nenberger, P. H.; McCammon, J. A. J. Mol.
Biol. 1999, 285, 1811.
(48) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger,
P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular
simulation: The GROMOS96 manual and user guide; Verlag der Fachvereine:
Zu¨rich, 1996.
(49) Daura, X.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem.
1998, 19, 535.
(50) Scott, W. R. P.; Hu¨nenberger, P. H.; Tironi, I. G.; Mark, A. E.;
Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Kru¨ger, T.; van
Gunsteren, W. F. J. Phys. Chem. A 1999, 103, 3596.
(51) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.;
Hermans, J. In Intermolecular Forces; Pullman, B., Ed., Reidel: Dordrecht,
1981; pp 331-342.
(52) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Haak,
J. R. J. Chem. Phys. 1984, 81, 3684.
(53) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys.
1977, 23, 327.
(54) Madura, J. D.; Davis, M. E.; Gilson, M. K.; Wade, R. C.; Luty, B.
A.; McCammon, J. A. In ReViews in computational chemistry; Lipkowitz,
K. B., Boyd, D. B., Eds.; VCH Publishers, Inc.: New York, 1994; pp 229-
267.
(55) Chakrabartty, A.; Baldwin, R. L. AdV. Protein Chem. 1995, 46,
141.
Molecular Dynamics Simulations of an Octapeptide J. Phys. Chem. B, Vol. 104, No. 15, 2000 3675