| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Biophys J, February 2000, p. 668-682, Vol. 78, No. 2
and
*Department of Biochemistry, Virginia Polytechnic Institute and
State University, Blacksburg, Virginia 24061, and
Laboratory of Structural Biology, National Institute of
Environmental Health Science, Research Triangle Park, North Carolina
27709 USA
 |
ABSTRACT |
|---|
 TOP
 ABSTRACT
 INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
Molecular dynamics (MD) simulations of the DNA duplex
d(CCAACGTTGG)2 were used to study the relationship between
DNA sequence and structure. Two crystal simulations were carried out;
one consisted of one unit cell containing two duplexes, and the other
of two unit cells containing four duplexes. Two solution simulations were also carried out, one starting from canonical B-DNA and the other
starting from the crystal structure. For many helicoidal parameters,
the results from the crystal and solution simulations were essentially
identical. However, for other parameters, in particular, 
, 
, 
,
(
), phase, and helical twist, differences between
crystal and solution simulations were apparent. Notably, during crystal
simulations, values of helical twist remained comparable to those in
the crystal structure, to include the sequence-dependent differences
among base steps, in which values ranged from 20° to 50° per base
step. However, in the solution simulations, not only did the average
values of helical twist decrease to ~30° per base step, but every
base step was ~30°, suggesting that the sequence-dependent
information may be lost. This study reveals that MD simulations of the
crystal environment complement solution simulations in validating the
applicability of MD to the analysis of DNA structure.
| |
INTRODUCTION |
|---|
|
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
Understanding factors associated with essential
cellular processes such as replication and transcription will lead to
an enhanced ability to regulate these processes and thereby treat a
variety of diseases more effectively. Of importance in replication and transcription are the structure and dynamics of DNA, which are dependent to some extent on DNA sequence (Hunter, 1996
). Thus an
improved understanding of the relationship between DNA sequence and
structure is a subject of considerable interest.
Methods that may be used to examine the relationship between DNA
sequence and structure include x-ray crystallography, NMR, and
theoretical methods such as molecular dynamics (MD). X-ray crystallography has proved to be particularly valuable in the examination of DNA structure, but the principal limitation of this
technique is the need to obtain high-quality crystals for every
sequence of interest. Another criticism that is sometimes invoked is
that crystal packing effects may constrain the DNA structure, thereby
confounding an analysis of the relationship of sequence to structure
(Ramakrishnan and Sundaralingam, 1993; Tippin and Sundaralingam, 1997).
NMR has also proved valuable in the structural analysis of DNA, but a
problem associated with NMR is that the data are relatively
short-ranged, and consequently fewer data can be collected when an
extended structure like DNA is studied, as compared to more globular
structures like proteins (Hartmann and Lavery, 1996). As a result, NMR
structures of DNA, in the absence of dipolar anisotropic coupling data,
tend to be of a lower resolution than crystal structures. Theoretical
approaches such as MD are being applied more frequently in studies of
DNA structure, especially since the application of Ewald methods to electrostatics has led to better agreement between simulation and
experiment (Darden et al., 1999; York et al., 1993, 1994, 1995).
However, MD is still not sufficiently tested to be used routinely to
study sequence-dependent structural differences in DNA. Nevertheless,
it is important to extend the applicability of MD because it can
provide an atomic-level description of DNA structure and dynamics that
will prove valuable for interpreting experimental data and developing
models of DNA structure and function.
Recent review articles highlight various aspects of MD simulations of
nucleic acids (Cheatham et al., 1998; Jayaram and Beveridge, 1996;
Westhof et al., 1995). Several recent papers describe results that are
pertinent to our studies and that exemplify the current state of MD
simulations of nucleic acids. For example, Cheatham and Kollman (1996)
performed MD on DNA duplexes of the sequence d(CCAACGTTGG)2, as we have done. They carried out
four MD simulations, two starting from canonical B-DNA and two from
canonical A-DNA. All simulations produced structures with
characteristics of B-DNA that were within 0.8-1.6 Å root mean square
(RMS) deviation of one another and 3.1-3.6 Å RMS deviation of the
published crystal structure (Prive et al., 1991). During the MD
simulation starting from canonical A-DNA, a transition to B-DNA was
observed on a nanosecond time scale, indicating that the energy barrier
for transition from A- to B-DNA was not prohibitive. Other recent studies include MD simulations of the dodecamer
d(CGCGAATTCGCG)2. This structure has been studied
extensively because it was the first duplex to be crystallized and
because it contains the recognition site of the EcoRI
restriction enzyme. In one of the studies, Duan et al. (1997) ran
simulations starting from the crystal structure of the duplex as it
exists in a complex with the endonuclease. In this complex the duplex
is kinked. MD simulation produced a structure that was no longer kinked
and was closer, based on RMS deviation, to the crystal structure of the
duplex in the absence of the protein. These investigators also noted
some evidence of sequence-dependent structural features and a spine of
hydration in the minor groove. Beveridge and co-workers have carried
out extensive MD simulations on nucleic acids, including two very recent studies on the d(CGCGAATTCGCG)2 dodecamer
(Young et al., 1997a,b). The starting structure for these studies was
canonical B-DNA. Several simulations were carried out to examine the
effects of initial counterion placement and the method of calculating electrostatic interactions (i.e., truncated potentials versus particle
mesh Ewald). One of the significant findings of these studies was that
Na+ counterions may enter the spine of hydration
in the minor groove of DNA, possibly in a sequence-specific manner.
This observation extends hypotheses that water molecules that
constitute the spine of hydration may stabilize DNA structure.
An impressive aspect of the simulations cited above is that structural
features characteristic of B-DNA are maintained over several
nanoseconds. However, not all features are preserved (e.g., helical
twist), so it is necessary to reconcile the results from simulation
with those from experimental methods such as crystallography. An
additional approach that can be applied in MD is simulations that start
from the crystal structure, with the crystal environment maintained
throughout the simulation (crystal simulation). Differences that arise
between crystal structure and simulated structure may facilitate the
identification of force-field parameters that need to be adjusted to
improve agreement between experiment and simulation. Some crystal
simulations of nucleic acids have been reported. For example, Darden
and co-workers performed a series of crystal simulations when they were
demonstrating the efficacy of the Ewald methods. An MD simulation of
the dodecamer d(CGCGAATTCGCG)2 in the crystal
unit cell yielded a stable trajectory and an average structure with an
RMS deviation for all heavy atoms of 1.16 Å from the crystal structure
(York et al., 1995). Similarly, an MD simulation of Z-DNA
(d[CGCGCG]2) in its crystal environment produced an average structure with an RMS deviation of 0.5 Å from the
crystal structure (Lee et al., 1995b). MD simulations of the RNA
dinucleotides ApU and GpC were also undertaken (Lee et al., 1995a).
These molecules represent a particular challenge for MD because all
heavy atom positions in the crystal are known, to include nucleotide,
counterion, and water positions. The average structures from MD
simulations of the dinucleotides were within 0.4 Å of the crystal
structures, and the calculated and experimental temperature factors
were comparable. These results suggest that crystal simulations may
provide a means of performing MD simulations under conditions in which
many structural features of nucleic acids are preserved without the
molecules being overly constrained within the crystal lattice.
Moreover, a direct comparison between crystal structure and crystal
simulation provides a way to rigorously evaluate the quality of the simulations.
As noted above, crystal packing effects are sometimes cited as a
limitation of crystallography in the analysis of DNA structure. Dickerson has effectively argued that crystal packing may constrain the
DNA conformation in the crystal, but the DNA can only adopt structures
that are favored by its sequence (Dickerson et al., 1994). That is,
crystal packing does not force the DNA into structures that it could
not otherwise assume. It also bears mentioning that analysis of
structures in solution, whether by NMR or MD, also may not represent a
realistic environment for the state of DNA in cells. In the nucleus,
where the density is 1.3-1.4 g/ml (Sheeler and Bianchi, 1980), DNA is
tightly packed and constrained in a manner than may be more similar to
its crystalline than to its solution state.
In the studies described here, we report MD simulations of the DNA decamer d(CCAACGTTGG)2 in the crystal environment. A simulation including one unit cell containing two duplexes was run for 25 ns, and a simulation including two unit cells containing four duplexes was run for 15 ns. We also conducted solution simulations for comparison. These long simulations enabled us to examine in detail the stability of the duplexes during the simulations.
| |
METHODS |
|---|
|
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
General methods
Four MD simulations were conducted as part of this project. Two of
them were simulations of the crystal environment, one involving a
single unit cell containing two DNA duplexes, and the other two unit
cells containing four duplexes. The other two simulations were solution
simulations in which a single duplex was immersed in a box of water.
For one of the solution simulations the starting structure was
canonical B-DNA, and for the other the starting structure was the
crystal structure. Structures were visualized using Visual Molecular
Dynamics (VMD) (Humphrey et al., 1996). All modeling was done with
AMBER4.1 programs (Pearlman et al., 1995) implemented on Silicon
Graphics computers with R10000 processors. The Cornell et al. (1995)
force field was used. MD simulations were conducted with the SANDER
module of AMBER. Simulations were carried out with a 2-fs time step at
300 K. The SHAKE algorithm was applied to all bonds involving hydrogen
atoms. Van der Waals interactions were calculated using an 8-Å
atom-based nonbond list, while the long-range Coulomb energy was
evaluated by the particle-mesh-Ewald (PME) method (Essmann et al.,
1995). Specific details pertaining to each simulation are given below.
Crystal simulations
One unit cell
The starting coordinates for the crystal simulations were taken from the PDB file with code 5dnb (Prive et al., 1991). The PDB file
contains coordinates for only one strand of the duplex, so the
remaining strand was generated through the appropriate symmetry
transformation. A second duplex was added to the unit cell by
translating the initial duplex according to guidelines published by
Prive et al. (1991). An illustration of the arrangement of the two
duplexes is given in Fig. 1. The
crystallographic positions of the Mg2+ ions were
maintained during the generation of the unit cell, with 14 Mg2+ ions present in the system. Crystallographic
waters were added by applying the appropriate symmetry transformations,
as was done to generate the duplex structures, and then by adding
hydrogen atoms, using the GWH utility within AMBER. The EDIT and PARM
modules of AMBER were then used to generate topology and coordinate
files containing positions of DNA atoms, Mg2+
ions, and crystallographic water molecules. Bulk water molecules were
then added to the system, and minimization was carried out to adjust
the positions of water molecules that were in close contact with other
atoms in the system. To achieve electroneutrality, eight bulk water
molecules were selected randomly and converted into
Na+ ions. The energy of the system was then
minimized to remove bad atom-atom contacts. After minimization, MD was
run for 100 ps, during which the DNA, Mg2+ ion,
and crystallographic water positions were fixed and the bulk water
molecules and Na+ ions were allowed to move. We
then attempted to perform MD with decreasing constraints on the
crystallographic waters, but we encountered unrealistically high
energies. This problem was traced to the wrapping option in the SANDER
module of AMBER. That is, water molecules that were wrapped to an
adjacent periodic box were being constrained to a position in the
original box, thereby generating large constraint energies. Thus
wrapping was turned off for subsequent MD simulations, which would have
no effect on forces and energies but was done for convenience of
analysis. MD was then carried out at 300 K, during which bulk water
molecules and Na+ ions were allowed to move, but
crystallographic waters were constrained so that they did not move a
significant distance from their crystallographic positions. The system
was then minimized, first by minimizing bulk waters and
Na+ ions; then bulk waters,
Na+ ions, and crystal waters; then all hydrogen
atoms; and finally the complete system.
|
Two unit cells
A system consisting of two unit cells was constructed by translating a copy of the original unit cell in the Z-dimension (Prive et al., 1991), giving a system of 5764 atoms. This simulation enabled
us to increase the amount of sampling during a single simulation and to
determine if the results were altered by relieving the end-to-end
constraints imposed by periodic boundary conditions. MD was run by
heating the system from 0 to 300 K over 30 ps, running constant-volume
MD for an additional 50 ps, and then running constant-pressure MD for a
total simulation time of 15 ns.
Solution simulations
Canonical B-DNA
This solution simulation was based on methods described by Cheatham and Kollman (1996). Canonical B-DNA was built using the NUCGEN
module of AMBER. Na+ counterions and water were
added in LEAP. The size of this system was 9440 atoms. The system was
equilibrated by running MD for 80 ps at 500 K and then for 20 ps at 300 K. The system was then subjected to energy minimization, first with the
DNA atoms fixed and then with no constraints. After minimization, the
system was heated first by running MD for 10 ps at 100 K, then for 10 ps at 200 K, and finally for 23 ps at 300 K, with all of these steps at
constant volume. Constant pressure conditions were then imposed for the
remainder of the simulation, which proceeded for a total of 12 ns.
Crystal structure
The crystal coordinates from PDB file 5dnb were used for the DNA structure. The crystallographic counterions and water molecules were not used. Na+ counterions and water were added in LEAP, giving a system of 9083 atoms. The system was equilibrated by running MD for 100 ps while the DNA atom positions were fixed. The system was then minimized, by minimizing water and Na+ ions and then all atoms. The system was heated from 0 to 300 K over a period of 30 ps under constant-volume conditions. It was further equilibrated for 50 ps at constant volume, after which constant pressure was imposed for the remainder of the simulation, which was run for 5 ns.Data analysis
Coordinates for analysis were saved every picosecond during the
simulations. DNA structural parameters were calculated with Curves 5.1 (Lavery and Sklenar, 1988). A global helical analysis was used in the
Curves analysis. Except where noted, structural parameters for the
various simulations were obtained by averaging over the following time
periods: one unit cell crystal simulation, 20-25 ns; two unit cell
crystal simulation, 10-15 ns; solution simulation starting from
canonical B-DNA, 10-12 ns; solution simulation starting from the
crystal structure, 4-5 ns. Average values for helicoidal parameters
were obtained by averaging the values for the individual structures
collected at 1-ps intervals over the time periods indicated.
| |
RESULTS |
|---|
|
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
Analysis of DNA structure
Many of the data pertaining to DNA structure from the simulations are summarized in Table 1. Results from the one unit cell simulation will be discussed first, followed by the results from the two unit cell simulation, after which a comparison of these crystal simulations with the solution simulations will be made.
|
One unit cell simulation
Shown in Fig. 2 is the RMS deviation, relative to the starting structure, for the two duplexes in the one unit cell simulation. This long simulation allows for a careful analysis of the stability of the trajectory. Notably, the RMS deviation of Duplex 1 was relatively stable during the last 7 ns of the simulation, although considerable fluctuation in RMS deviation was noted earlier in the trajectory. Duplex 2 was relatively stable for the last 15 ns of the simulation. The RMS deviations of the duplexes were very similar to one another by the end of the trajectory.
|
|
Two unit cell simulation
In Fig. 3 are the RMS deviations for the simulation involving four duplexes in two unit cells. During the first 2.5 ns, the RMS deviations were similar for all of the duplexes, and the system appeared to be relatively stable. Shortly thereafter, larger differences in RMS deviation were noted among the duplexes, and these differences persisted, even out to 15 ns. Examination of the helicoidal parameters for these duplexes (Table 1) reveals some differences, especially in the backbone torsion angles.
|
Comparison of crystal and solution simulations
The RMS deviations of the solution simulations suggested that a
stable trajectory was attained within 4 ns, which is faster than in the
crystal simulations (data not shown). It is possible that the crystal
environment imposes constraints that increase the time required to
achieve stability. The RMS deviations of the DNA duplexes in the
crystal and solution simulations relative to the crystal structure are
given in Table 2. It is seen that the structures from the solution
simulations deviated considerably more from the crystal structure than
did those from the crystal simulations. Notably, the values of RMS
deviation in the solution simulations are in the range of the RMS
deviations reported by Cheatham and Kollman (1996).
A more critical comparison of the structures can be made by considering
the helicoidal parameters listed in Table 1. The results for several of
the helicoidal parameters are particularly notable when the parameters
from the crystal and solution simulations are compared with one another
and with parameters for the crystal structure. These parameters include
, , , , , phase, and twist. Each of these parameters
will be examined in more detail.
and
from
the crystal simulations are consistently much lower than the value of
from the solution simulations or the crystal structure. These low
values of appear to be compensated for by the consistently higher
values of the torsion angle. The correlation between and is
more apparent from Fig. 4, in which
values of and from the solution simulation starting from the
canonical structure and the one unit cell crystal simulation are shown.
A number of points can be made regarding the data in Fig. 4. First,
changes in and appear to correlate well
as decreases, increases. Second, discrete changes in values of and are
observed during the simulations, suggesting that changes among distinct
structural forms are occurring. If a steady decline in and a steady
increase in had been observed over the course of the simulations,
we would have attributed the change to a force-field-related problem. Third, the differences in values reported in Table 1 are dependent on
the time frame used to calculate the average structural properties. For
the solution simulation starting from the canonical structure, we used
10-12 ns; for the one unit cell crystal simulation, we used 20-25 ns.
Values of during the crystal simulation were lowest during that
time period. During the 10-12-ns time period in the solution
simulation, values of were relatively high, and lower values were
observed during the 5-8-ns time period.
|
and phase angle
and the pseudorotation phase angle, the crystal
simulations gave much better agreement with the crystal structure than did the solution simulations. It is generally accepted that the torsion, being part of the deoxyribose ring, correlates with the
pseudorotation phase angle (Dickerson, 1992). Thus we will not discuss
in our more detailed analysis but rather consider the
pseudorotation phase angle, which includes as one of its components. We analyzed changes in the pseudorotation phase angle, which provides an analysis of the sugar repuckering, over the course of
the trajectories. From Fig. 5 it is
seen that during the simulations, the sugar conformation is most
frequently C2'-endo (phase angle = 137°-194°),
although transitions to C3'-endo (phase angle = 1°-34°) are observed. Repuckering tends to occur more often
during the solution simulation than in the crystal simulation. In fact,
the example chosen for illustration of the crystal simulations showed
the greatest frequency of transition among the duplexes in either of
the crystal simulations, although some repuckering was observed in all
of the duplexes. As is apparent from Fig. 5, when transitions occurred
during the crystal simulations, the sugar rings tended to stay in the
newly adopted conformation for an extended period, unlike transitions
in the solution simulations. This observation of less frequent
repuckering during the crystal simulation would explain the close
agreement between the value of the phase angle in the crystal
simulations and the crystal structure. That is, if frequent repuckering
to C3'-endo were to occur, as in the solution simulation,
then the average value of the phase angle would be correspondingly
lower. Moreover, the value of the phase angle reported in Table 1 is an
average over all residues. In the crystal simulation, the repuckering
occurred principally in one residue, T-8, and the phase for other
residues in the crystal simulations tended to be higher than that for
the corresponding residues in the solution simulation. These
differences in repuckering between the solution and crystal simulations
may reflect constraints imposed by crystal packing, although it is clear that repuckering can still occur during the crystal simulations. These crystal packing constraints may overcome force-field deficiencies such that better agreement with the crystal structure is obtained.
|
( )
for the structures in
the crystal simulations are generally somewhat higher than that for the
crystal structure, and the values in the solution simulations are
comparably lower. Values of in the crystal simulation are similar
to that in the crystal structure, but those in the solution simulations
are considerably higher. The difference ( ) is often used
when these values are analyzed, because it assumes values that are
characteristic of the BI backbone conformation (90°) and the BII
backbone conformation (+90°) (Lavery, 1994). In the crystal
structure, the values of ( ) for C-2 and T-8 are
indicative of the BII conformation, while the others are characteristic of the more common BI conformation. Notably, in the solution
simulations, all values are indicative of the BI conformation. In the
crystal simulations, some of the duplexes retain values of ( ) at positions C-2 and T-8 that are characteristic of the BII
conformation (data not shown).
We have examined values of ( ) as a function of time.
Shown in Fig. 6 are results from the
solution simulation starting from the canonical B-DNA structure, and in
Fig. 7 are results from one of the
crystal simulations. The results from the solution simulation
illustrate that transitions between BI and BII occur, albeit relatively
infrequently. This is most apparent at residue A-4. In the crystal
simulation, frequent transitions are observed at several bases,
including C-2 and T-8. Note that the example in Fig. 7 was chosen
because it illustrates frequent transitions. For other DNA strands in
the crystal simulations, the transitions were not as frequent. However,
it is important to recognize that the crystal environment does not
overly constrain the DNA duplexes such that they remain rigid in the
crystal lattice during the simulation. Rather, flexibility is apparent,
based on the transitions illustrated here.
|
|
Helical twist
Perhaps most interesting among the helicoidal parameters are the values for helical twist. Solution simulations of DNA with the Cornell et al. (1995) force field typically result in an undertwisting of
the duplex, such that values for twist are lower than those observed in
crystal structures (Cheatham and Kollman, 1996, 1997; Cieplak et al.,
1997; Young et al., 1997b). This phenomenon also is apparent in
our results (Table 1). In contrast, the twist parameter values for the
crystal simulations are higher than for the solution simulation and are
very close to that in the crystal structure. This result may derive
from the constraints placed upon the DNA duplexes within the unit cell.
Helical twist is also of interest because it appears to be sequence
dependent (El Hassan and Calladine, 1997; Gorin et al., 1995). For some
of the other parameters that are reported to show some sequence
dependence (e.g., roll, slide, propeller), the standard deviation
(fluctuation) over the course of the simulations was sufficiently large
that we could not discern a sequence dependence. Moreover, some of the
sequence dependence of DNA structural features deduced from crystal
structures has recently been attributed to experimental error in the
crystal structure determinations (Shui et al., 1998a). In our
simulations, the values of helical twist are sufficiently large and the
fluctuations sufficiently small that some sequence dependence may be
apparent. To examine the sequence dependence of helical twist in the
simulations, values were determined for each base step (Table
3). Notably, the values of twist at
each base step in the crystal simulations correspond well to the values
in the crystal structure. That is, the sequence dependence of helical
twist was maintained in the crystal simulations. However, values of
helical twist in the solution simulations were uniformly close to
30°. Because the sequence dependence of helical twist was not
preserved in the solution simulations, we examined the time course over
which the helical twist changed in the solution simulation starting
from the crystal structure. In Fig. 8 is
shown the decrease in helical twist for the three base steps that have the highest values in the crystal structure. The helical twist for two
of the base steps decreased to around 30° within the first 250 ps of
the trajectory, and the helical twist for the third base step decreased
within 1 ns.
|
|
Minor groove width
Another parameter that is reported to exhibit some sequence dependence is minor groove width (Boutonnet et al., 1993; Nelson et
al., 1987; Sarma et al., 1997). A comparison of minor groove widths
from the simulations with that of the crystal structure is given in
Fig. 9. The minor groove width is not
precisely preserved in any of the simulations. However, it is apparent
that a narrowing of the minor groove as observed in the crystal
structure is seen in the crystal simulations but not the solution
simulations. The time dependence over which the minor groove width
increased in the solution simulation starting from the crystal
structure was determined (Fig. 10).
Considerable fluctuation in minor groove width was observed during the
trajectory, but it appears to be stabilizing in its wider state by ~4
ns into the simulation.
|
|
Analysis of B-factors
Another means of comparing the results of MD simulation with the crystal structure is to calculate B-factors and compare them to the crystallographic B-factors. The B-factors provide a measure of the fluctuation at atomic positions within the DNA molecules. In Fig. 11, calculated B-factors from the one unit cell simulation and the solution simulation starting from the canonical structure are compared to the B-factors from crystallography. For clarity, each strand of DNA from the one unit cell simulation is depicted in a separate panel (Fig. 11, A-D). As expected, the regions of highest atomic fluctuation are in the backbone atoms. Agreement between calculated and experimental B-factors varies among the four strands in the one unit cell simulation. The closest agreement is observed in Fig. 11, B and D. In Fig. 11, A and C, the calculated results suggest regions of less atomic fluctuation than in the experiment. We also calculated B-factors from a solution simulation (Fig. 11 E). Clearly, the magnitude of the fluctuations is greater in the solution simulation than in the crystal simulation or in the experiment. However, the agreement between regions undergoing greater and lesser atomic fluctuation is good when calculated and experimental values are compared.
|
Mg2+, water structure, and dynamics
It also was of interest to consider the structural relationships
among the Mg2+ ions, water molecules, and DNA
during the crystal simulations to determine if these components behaved
as expected. The water molecules that constituted the hydration sphere
around each Mg2+ ion were examined for the
crystal structure and for structures at various times during the one
unit cell crystal simulation. In only two cases was an exchange of
water molecules observed within the primary coordination sphere. That
is, the water molecules constituting the primary sphere of hydration
tended not to exchange with the bulk solvent. The radial distance
distribution of water molecules relative to Mg2+
ions was determined for the last 5 ns of the one unit cell trajectory (Fig. 12). It is clear that the
Mg2+ ions maintain the hydration shell, with the
molecules in the inner sphere 2.0-2.1 Å from the
Mg2+ and those in the outer sphere 4.1-4.3 Å from the Mg2+ (Bock et al., 1994).
|
Analysis of the trajectories revealed that some of the Mg2+ ions (with their associated sphere of hydration) moved relative to the DNA molecules. An example of the extent of this movement is given in Fig. 13, in which the position of a Mg2+ ion relative to nearby DNA phosphorus atoms is illustrated. The motion is characterized by small fluctuations around equilibrium positions with a few larger abrupt changes in position. When viewed as a radial distribution function, the Mg2+-P distances show a primary distance of ~4.9 Å, which is characteristic of a system in which the Mg2+ hydration sphere is maintained and no direct Mg2+-P bonding occurs (data not shown).
|
| |
DISCUSSION |
|---|
|
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
Our crystal simulations appear to provide a reasonable
representation of the crystal environment in the unit cell. Under
constant-pressure MD conditions, the volume of the unit cell during the
last 5 ns of the one unit cell simulation was 26,286 ± 140 Å3, which agrees reasonably well with the
published value of 25,980 Å3 (Prive et al.,
1991). That the box size stayed essentially constant during MD
indicates that we had included the appropriate number of bulk water
molecules, in addition to the DNA, Mg2+
ions, and crystallographic water molecules, to simulate the
crystal environment.
One of our simulations, the solution simulation starting from the
canonical B-DNA structure, replicates the work of Cheatham and Kollman
(1996), though our simulation was carried out for a considerably longer
period of time (12,000 ps versus 1456 ps). Our results for mean values
of helicoidal parameters are in excellent agreement with those of
Cheatham and Kollman (1996), although our standard deviations typically
were considerably lower. Although we averaged over a longer time period
(2000 ps versus 856 ps), it is unlikely that the magnitude of this
difference is large enough to explain the different values for standard
deviations. The difference in standard deviations most likely results
from different methods of data analysis. For example, in our
calculations reported in Table 1, we calculated an average value for
each parameter in each structure of the trajectory and then calculated a mean and standard deviation based on these average values from all of
the structures. An alternative is to calculate averages of individual
residues, base pairs, or base-pair steps (Cheatham and Kollman, 1996).
Nevertheless, a careful comparison of our data with those of other
investigators indicates that the range of values observed over the
course of a trajectory is similar. Thus, we conclude that our solution
simulation conditions and results are comparable to those reported previously.
When comparing the helicoidal parameters from our crystal and solution
simulations with one another and with the crystal structure, we see
that differences are most apparent in the parameters , , ,
, phase, and helical twist. For many of the other parameters, the
absolute values are close to zero and the standard deviations are
sufficiently high that differences cannot be discerned. The correlated
crankshaft motion involving and is expected in that it has been
observed in other MD simulations (Swaminathan et al., 1991; Cheatham
and Kollman, 1997) and in NMR studies (Xu et al., 1998). As noted
above, the abrupt changes in and suggest that rapid transitions
among similar structures are occurring and that all structures are
energetically favorable. Notably, significant changes in RMS deviation
were not observed when and changed abruptly, suggesting that in
terms of the overall structure, the small changes in were balanced
by small changes in . Differences between values of and in
the crystal simulations and those in the solution simulations would
suggest some constraint due to crystal packing forces, while
differences in these values relative to the crystal structure would
suggest some deficiency in the force field. This problem with the force
field is also apparent from the observation that the lower average
values of and higher values of relative to the crystal
structure reflect to some extent the time period over which structures
were averaged. From Fig. 4 we can see that if averaging in the one unit
cell crystal simulation were to be carried out over the 10-12-ns time period, as was done for the solution simulation, values of and would be much closer to those from the solution simulation and the
crystal structure.
From analysis of the time dependence of and , as well as other
parameters, we do not think that the DNA was overly constrained during
the crystal simulations. For example, BI to BII backbone transitions
were apparent in the crystal simulations, just as in solution
simulations. Sugar repuckering also occurs in the crystal simulation,
although the rings near the ends of the duplex appear to be more
tightly constrained than those toward the middle of the duplex, a
finding that is not observed in the solution simulations. From the
B-factors, it also appears that DNA in the crystal simulation is
somewhat less flexible than that in the solution simulations,
suggesting that the tighter packing within the simulated crystal is
constraining movement to some extent. However, the close agreement
between B-factors from the crystal simulations and those from the
crystal structure suggests that the conditions for our crystal
simulations, along with the force field of Cornell et al. (1995),
accurately represent the dynamics of the duplex in the crystal environment.
Our analysis of helical twist and minor groove width make apparent some
problems associated with solution simulations of DNA. Other
investigators also have noted that the average value for overall
helical twist decreases to ~30° during solution simulations (Cheatham and Kollman, 1996, 1997; Cieplak et al., 1997; Young et al.,
1997b). The dependence of helical twist on base step in MD simulations
has been less thoroughly analyzed. Duan et al. (1997), in a 1-ns
simulation of the dodecamer duplex
d(CGCGAATTCGCG)2, reported an unwinding of the
duplex and a convergence of twist angles to ~33° within the
first 500 ps of their 1-ns simulation. The average value of
helical twist may decrease even further in a longer simulation. In any
case, our result is consistent with theirs in that differences in
helical twist as a function of base step are not apparent in solution
simulations. The loss of this sequence-dependent information is a
problem that must be resolved if solution MD simulations are to be used
in a predictive way to understand the relationship between DNA sequence
and structure. That we were able to preserve this sequence-dependent
information, even in extremely long crystal MD simulations, suggests
that these crystal simulations may provide a tool for analyzing the
variables within the solution simulations that account for the decrease in helical twist. Moreover, Cheatham et al. (1999) have recently released a modified version of the force field of Cornell et al. in
which improved values for pseudorotation phase angle and helical twist
are observed in solution simulations. We intend to apply this modified
force field to our simulations in the future.
The minor groove width also shows sequence dependence (Boutonnet et
al., 1993; Nelson et al., 1987; Sarma et al., 1997) and thus merits
careful scrutiny in MD simulations. We noted that some narrowing of the
minor groove is apparent in all of the duplexes in our crystal
simulations, although we clearly have not maintained the minor groove
width of the crystal structure precisely. In our solution simulations,
little narrowing of the minor groove is observed. In the simulations of
Cheatham and Kollman (1996, 1997) of the
d(CCAACGTTGG)2 decamer, the minor groove width
was narrowed, but not to the same extent as in the crystal structure. In the study by Duan et al. (1997) of the Dickerson dodecamer starting
from the "bent" crystal structure, the duplex appeared to retain
some narrowing within the minor groove, although the simulation was
performed for only 1 ns. In a 5-ns simulation of the Dickerson
dodecamer starting from canonical B-DNA, Young et al. (1997a) noted a
slight narrowing of the minor groove, though not to the same extent as
seen in the crystal structure. The differences observed among these
studies probably result from the different sequences under
investigation and from different simulation conditions. In particular,
as we have illustrated in Fig. 10, the minor groove width changes with
time, with a narrow minor groove more likely to be observed on shorter
time scales (e.g., <1.5 ns) when we start from a crystal structure
with a narrow minor groove.
The behavior of the Mg2+ ions during the crystal
simulations also needs to be considered because it is recognized that
metal ions may have a significant influence on DNA structure (Saenger, 1984). In our crystal simulations, the Mg2+ ions
from the crystal structure were included, with eight
Na+ ions needed to achieve electroneutrality. We
were most interested in the behavior of the Mg2+
ions during the simulations and analyzed their positions relative to
the DNA as a function of time. In addition, Na+
ions generally are considered to be very diffusible, although they are
reported to form inner sphere complexes with DNA (Westhof et al.,
1995). Recent studies support this observation. For example, Young et
al. (1997a) observed that Na+ ions are
incorporated within the primary hydration sphere of the minor groove of
d(CGCGAATTCGCG)2, from which they suggest that
Na+ ions may function as an extrinsic source of
sequence-dependent effects on DNA structure. This observation was
reinforced by studies of Shui et al. (1998a,b). In a high-resolution
crystal structure of d(CGCGAATTCGCG)2, they noted
that Na+ ions penetrated the spine of hydration
in the minor groove (Shui et al., 1998a). They subsequently confirmed
this observation in a structure of the dodecamer in the presence of
K+ ions, which are easier to distinguish from
water molecules than are Na+ ions (Shui et al.,
1998b).
The coordination of divalent cations is hypothesized to occur in two
ways. First, an outer sphere complex can form; then, with the loss of
one or more water molecules, an inner sphere complex can form (Westhof
et al., 1995). However, Mg2+ ions are reported to
form inner sphere complexes at a relatively low rate
(~105 s
1) and primarily
with short polyadenylate chains (Porschke, 1978). Cowan et al. (1993),
using 25Mg NMR spectroscopy, determined that
binding of Mg2+ ions to DNA involves outer sphere
complexes with
Mg(H2O)62+
and that binding to G/C-rich sequences occurs with 40- to 100-fold higher affinity than binding to A/T-rich sequences. Similarly, Shui et
al. (1998b) made a survey of complexes between
Mg2+ ions and B-DNA in the Nucleic Acid Data Bank
(Berman et al., 1992) and noted that when Mg2+
ions are bound to DNA, they remain fully hydrated and interact with DNA
through water molecules in the first hydration shell. However, in a
recent high-resolution (1.1 Å) crystal structure of
(CGCGAATTCGCG)2 most of the contacts between
Mg2+ ions and DNA were mediated by water
molecules, but two direct contacts between the ions and DNA were
observed (Tereshko et al., 1999). Moreover, when Buckin et al. (1994)
applied ultrasonic titration experiments in the analysis of
Mg2+-DNA binding, they concluded that
Mg2+ ions form outer sphere complexes with
oligomers containing dA-dT base pairs and inner sphere complexes with
dG-dC base pairs.
In MD simulations in which Mg2+ ions are present,
the ions exist primarily as hexahydrates (Lee et al., 1995b;
MacKerrell, 1997). Considering the length of the simulations and the
relatively slow exchange rate for H2O bound to
Mg2+ ions noted above, it is likely that inner
sphere coordination would not be observed. However, in the simulations
of Lee et al. (1995b), two of the Mg2+ ions were
hypothesized to exist in a complex in which the ions were bridged by
hydroxide ions. Clearly, additional studies are required to understand
binding of Mg2+ ions to DNA and the effect on DNA
structure, which is a research problem that we continue to pursue.
| |
ACKNOWLEDGMENTS |
|---|
DRB thanks the National Institute of Environmental Health Science (NIEHS) for partial support of a sabbatical year. LGP thanks the North Carolina Supercomputing Center for a grant of computer time and the National Institutes of Health for grant HL-06350. This research was performed while DRB was on study-research leave at the NIEHS.
| |
FOOTNOTES |
|---|
Received for publication 7 June 1999 and in final form 8 October 1999.
Address reprint requests to Dr. David R. Bevan, Department of Biochemistry, 201 Fralin Biotechnology Center, Virginia Tech, Blacksburg, VA 24061-0308. Tel: 540-231-5040; Fax: 540-231-9070; E-mail: drbevan{at}vt.edu.
| |
REFERENCES |
|---|
|
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES
|
|---|
visual molecular dynamics.
J. Mol. Graph.
14:33-38[Medline].ApU and GpC.
Chem. Phys. Lett.
243:229-235.a comparative study of the isoforms of the tetragonal and hexagonal family of octamers with differing base sequences.
J. Biomol. Struct. Dyn.
11:11-26[Medline]. and with 13C chemical shifts.
J. Am. Chem. Soc.
120:4230-4231
Biophys J, February 2000, p. 668-682, Vol. 78, No. 2
© 2000 by the Biophysical Society 0006-3495/00/02/668/15 $2.00
This article has been cited by other articles:
 |
|
 |
|
  P. K. Maiti, T. A. Pascal, N. Vaidehi, J. Heo, and W. A. Goddard III Atomic-Level Simulations of Seeman DNA Nanostructures: The Paranemic Crossover in Salt Solution Biophys. J., March 1, 2006; 90(5): 1463 - 1479. [Abstract] [Full Text] [PDF]
|
|
