J. Phys. Chem. 1995,99, 11591- 11599
11591
Counterion Distribution around DNA Studied by Molecular Dynamics and Quantum Mechanical Simulations C. A. Laughton The CRC Biomolecular Structure Unit at The Institute of Cancer Research, Sutton, Surrey SM2 5NG, England, UK
F. J. Luque Departament de Farmacia, Unitat Fisicoquimica, Facultat de Farmacia, Universitat de Barcelona, Avgda Diagonal sn, Barcelona 08028, Spain
M. Orozco" Departament de Bioquimica, Facultat de Quimica, Universitat de Barcelona, C/ Marti i Franquks I , Barcelona 08028, Spain Received: October 31, 1994; In Final Form: March 21, 1995@
Quantum mechanical calculations and molecular dynamics (MD) simulations have been used to obtain an accurate and time-averaged representation of the charge distribution of a DNA-water-Na' system. Having derived a set of time-averaged charges for the phosphate and the sodium ions, two 150 ps MD simulations have been performed to obtain a more complete picture of the dynamic characteristics of the DNA, including solvent and counterion environment. In agreement with other theoretical approaches to the analysis of DNA/ ion interaction, about 80% of the counterions are associated with the DNA. This occurs almost entirely through interaction with the phosphate groups, the counterions being distributed approximately equally between the first and second coordination shells of the phosphates. The revised partial charges for the phosphate group and sodium ion appear to result in structural parameters for the DNA which remain in somewhat better agreement with experimental data than when the standard charges are used.
Introduction Molecular modeling, and particularly molecular dynamics (MD), has become a very powerful tool to complement experimental techniques for the study of DNA structure and functionality. The technique allows one to examine dynamical aspects of DNA structure inaccessible from X-ray crystallography and relate NMR data to structure with atomic detail.'-* The simulation of DNA structures is especially challenging, not only because of the size of the systems but also because of their chemical nature. At first sight, DNA structures are built of fragments which would not be expected to join together readily. The most obvious example of potential instability in DNA is the repulsion between phosphate groups. These groups, each with a net charge of - 1 e, are typically spaced every 7 8, along an oligonucleotide strand. The interaction of each with its two nearest neighbors implies a repulsive energy of almost 50 kcal/mol in the gas phase. Obviously, this repulsion, which would destroy the structure of DNA, is moderated in physiological media by the presence of water and, most importantly, positively charged counterions. The interaction between DNA and its counterion environment has been the subject of a large number of studies.'-20 The reasons for this interest are related not only to the importance of the phosphate-counterion interaction for the stability of DNA but also to the fact that it is a prototypical process whose rationalization would contribute to our understanding of the interactions between drugs and the DNA. Unfortunately, the study of counterion-DNA interactions is difficult from an experimental point of view due to the large mobility of @
Abstract published in Advance ACS Absrracts, June 15, 1995.
counterions. Thus, most of the studies have been done at the theoretical level; of particular note are the studies done with rigid representations of the DNA and Monte Carlo samplings of the counterions m o ~ e m e n t . ~ - ~Some ~ % ' ~studies have been done using explicit representations of ~ a t e r ? . ' ~but , ' ~most have been performed using a continuum representation of the aqueous en~ironment."-'~Very few studies have been performed using MD techniques on a fully flexible DNA system and an atomic model for Molecular dynamics, like Monte Carlo and molecular mechanics, is generally based on the use of a nonpolarizable classical Hamiltonian, the force field. The selection of the force field is probably the most critical aspect of any MD study and especially for a complex molecule like DNA. Studies by K o h a n and co-workers suggest that the AMBER all-atom force field (AMBEWAA)22can correctly reproduce the vibrational s p e ~ t r aof~ nucleic ~ , ~ ~ acid components, as well as the stacking and H-bond interaction of nucleic bases in both vacuum and aqueous s o l ~ t i o n . This ~ ~ , explains ~~ the wide use of this force field in the study of DNA structures and interactions. In AMBEWAA the atomic charges (key parameters for DNA simulations) were determined in vacuum using Momany's strategy,26typically at the STO-3G 27 leve1.22,28Considering the goodness of the correlation between STO-3G and 6-31G* charges and dipole^,^^-^' and the fact that the bases are buried inside the DNA structure, the vacuum STO-3G charges are expected to reproduce with reasonable accuracy the electrostatic properties of the nucleic bases and the deoxyribose. Unfortunately, this quality cannot be expected for the charges at the phosphate group, since (i) the basis set used to determine the charges at the phosphate group (STO-3G*) in AMBEWAA is
0022-365419512099-11591$09.0010 0 1995 American Chemical Society
Laughton et al.
11592 J. Phys. Chem., Vol. 99, No. 29, 1995
probably too small to reproduce the properties of a large anion (ii) the phosphate group in the DNA is embedded in water, in an environment where the charges in the gas phase may be no longer (iii) the counterion environment of DNA is expected to produce a notable polarization effect on the phosphate charge distribution, which is ignored in the AMBER/ AA charges; and finally (iv) the 3-D structure of DNA may introduce important asymmetries in the charges at the oxygens of the phosphate group, which are also neglected in the standard force field. In this paper, we describe a MD-based investigation of the interactions between DNA and Na+ counterions. For this purpose, previous to the MD simulations we performed a careful reparametrization of the phosphate-Na+ system to help overcome the problems noted above. The results of extended MD simulations with the revised parameters were compared with another conducted using the standard AMBEWAA force field. In order that these studies could be readily compared to previous work, they were performed on the fully hydrated DNA dodecamer d(CGCGAA'ITCGCG), for which extensive MD simulations have been r e p ~ r t e d ~ + and ' ~ +for ' ~ ?which ~ ~ a highresolution crystal structure is available.34 Methods Molecular Dynamics. Starting coordinates for the DNA heavy atoms were taken from the crystal structure,34 and hydrogens were added according to geometrical considerations. The system was made electrically neutral through the addition of sodium counterions; one being placed 4.0 8, from each phosphorus atom along the line bisecting the associated P-OA and P-OB bond vectors.35 Finally, the solute was placed in a box of Monte Carlo TIP3P water,36whose dimensions were such that no atom of the solute was closer than 8.0 A from any edge of the box. The final system consists of 24 nucleotides, 22 Na' ions, and 2675 water molecules, making a total of 8805 atoms. All MD simulations were performed under constantpressure periodic boundary conditions, and the simulation temperature (T = 298 K) was controlled through coupling to a heat bath.37 A time step of 0.002 ps was used in conjunction with and snapshots were saved every 0.2 ps. A 9 A nonbonded residue-based cutoff was used and the pair list updated every 20 dynamics steps, as standard in the AMBER package. Before allowing the DNA to move, some equilibration of the water and counterion environment was performed, and before allowing the counterions to move, the water structure was allowed to equilibrate. Therefore, for the first 20 ps of the dynamics the DNA and counterions were restrained to their initial positions with a force constant of 100 kcaV(mol/A2),and the temperature was held at 10 K for the first 3 ps and then raised to 300 K for the remaining 17 ps. For the next 20 ps only the DNA was constrained, permitting the counterions to partially equilibrate. Then five 2 ps dynamics runs were performed with the restraints on the DNA reduced from 100 kcal/(mol A) to 50,25, 10,5, and finally 1 kcal/(mol A). Only then was unrestrained dynamics on the whole system begun. Charge Parametrization. The determination of the phosphate-Na+ charge distribution in solvated B-DNA structures presupposes knowledge of the spatial arrangement of DNA, H20, and Na+, the determination of which in turn requires knowledge of the phosphate-Na+ charge distribution. Hence, an iterative self-consistent process was necessary (see Figure l), where MD and charge parametrization were repeated until convergence in the P-Na+ distribution. A first MD simulation (simulation 1) was performed with the standard AMBEWAA charges for 50 ps after equilibration
Initial charges
MDsimulation
Calculation of
-+ population factors
New populationweighted partial charges
Final charges
Figure 1. Flow diagram illustrating the iterative procedure for obtaining partial charges for the phosphate group and sodium counterion.
OAdP\ OB
Figure 2. Definition and illustration of the boundaries between the different families of phosphatekounterion configurations. TABLE 1: Definition of PhosphatdCounterion Conformational Families family la lb 2a 2b 3a 3b 4a 4b 5a 5b 6a 6b
7
description
within 4 8, of a P atom and within 30" of P - O m - O B bisector between 4 and 6 8, of a P atom and within 30" of P-ON P-OB bisector within 4 8, of a P atom and within 30" of P-OA vector between 4 and 6 8, of a P atom and within 30" of P-OA vector within 4 8, of a P atom and within 30" of P-OB vector between 4 and 6 8, of a P atom and within 30" of P-OB vector within 4 8, of a P atom and within 30"-50" of P-OA vector between 4 and 6 8, of a P atom and within 30"-50" of P-OA vector within 4 8, of a P atom and within 30"-50" of P-OB vector between 4 and 6 8, of a P atom and within 30"-50" of P-OB vector within 4 8, of a P atom but not in any of the above families between 4 and 6 8, of a P atom but not in any of the above families further than 6 8, from the nearest P atom
of the dodecamer-H20-Na+ system (see previous section). Analysis of the spatial arrangement of the phosphates and Na+ ions led to the identification of different regions of (phosphatedefined) space were the Na+ were most likely to be located. Then, selecting the P-Na+ distance and angles between the P-Na+ and P-OM-OB vectors and their bisector as descriptive variables (Figure 2 and Table l), we defined 13 conformational families, and the relative population of each was computed. Around 5-6 snapshots representing each populated family (11 of the 13 possible) were selected from the MD simulation to compute ab initio electrostatic charges.26 Due to the size of the system, a reduced model (the phosphate capped by two
Counterion Distribution around DNA methyls, the Na+, and the waters within 7 A) was used in all the quantum mechanical calculations. (This system is large enough to represent both ion-ion and ion-water-ion structures.) The methyl groups were placed in the positions corresponding to the C5' and C3' ribose atoms, and the 0 - C bond was rotated in such a way that the methyl hydrogens mimicked the ribose skeleton. The RHF wave functions of these systems were computed using a mixed basis set which would be expected to reproduce the charges of the system quite a ~ c u r a t e l yP: ~and ~ 0 were described at the 6-31+G(d) leve140 and the rest with the 6-31G(d) basis set4' All the waters were explicitly introduced into the RHF calculations as classical particles (TIP3P potential^^^), which interact with the phosphateNa+ system via a perturbational operator (see eq 1). M
3n
where 9 stands for the phosphate-Na+ gas phase Hamiltonian, n is a water molecule, and Qi are the charges at each of the three centers of the T P 3 P water model. The wave functions for each configuration of the phosphateNa+-water system were used to determine the electrostatic charges following a modified ~ e r s i o n * of * ~Momany's ~~ strategy.26 Thus, the quantum mechanical and classical molecular electrostatic potentials were fitted in two layers placed at 1.4 and 1.8 times the atomic van der Waals radii; a density of 5 points/A2 was used to guarantee the statistical quality of the charge^.^*-^' Wave funcions were determined using the HONDO-7 computer Molecular electrostatic potential and electrostatic charges were determined using MOPETEMOPFIT computer programs.43 The charges for every family were determined by averaging the charges for each representative configuration. The global charges for the phosphate-Na+ system were determined by a population-weighted averaging of the charges of the different families (see eq 2).
where j stands for each phosphate-Na+-water configuration, Pk stands for the population of family k during the dynamics, and {Qk} are the average charges of family k. Once the averaged charges were determined, a new MD simulation was performed (simulation 2), and after 50 ps the population analysis was repeated. This time, all 13 possible families were observed (see Table 1). Having calculated partial charges for the two extra families of configurations as described above and recalculated the population-weighted charges, the new averaged charges were used in a further 60 ps of MD simulation. The analysis of this new trajectory led to slightly altered population factors, and the averaged charges calculated using eq 2 also changed slightly. An additional 30 ps of MD did not introduce any major changes in the population factors, so the MD run with these third-iteration charges was extended to 150 ps (simulation 3A). To increase the amount of data for analysis and check to what extent the results were trajectory-dependent, a second 150 ps simulation (simulation 3B) was performed with this set of charges, starting from the same conformation as simulation 3A, but with altered velocities (see below). The population data were also used to compute the distribution functions for the phosphate-Na+ systems. This can be done using a generalized expression for the distribution function (eq 3), where the index i stands for a family of
J. Phys. Chem., Vol. 99, No. 29, 1995 11593 configurations, r is the distance P-Na+, and 8 is the angle between the line P-Na+ and the vector used to define that family (P-OA, P-OB, P-OAR-OB bisector; see Table 1 and Figure 2).
Both population and distribution functions can be used to obtain estimates of the free energy change due to changes in the position of Na+ ions with respect to the phosphate group. The free energy change due to the movement of a Na+ ion from one family to another at a constant P-Na+ distance can be obtained by integration of the distribution functions over the range of 8 corresponding to each family (eq 4). Furthermore, the total change in free energy due to the migration of the sodium ion from one family to another (changing both distance and orientation) can be obtained by integration over the range of 8 and r of the corresponding distribution functions (i.e., from the population numbers; see eq 5).
AGj+ = - R T l n ( h , B'+dB' gr(O')/he+deg;(e))
(4)
Free energy evaluations using eqs 4 and 5 are subject to important uncertainties due to problems in the equilibration of the phosphate Na+ distribution. To overcome or at least reduce this problem, we have followed the movement of each of the 22 phosphate-Na+ systems in two divergent (see below) MD simulations, each of 150 ps. Then, by the averaging of the information obtained for the 22 subsystems in the two 150 ps MD simulations, we have information in many ways equivalent to that from 6.6 ns of MD on a hydrated DNA/monophosphate/ Na+ system.@ Results and Discussion Simulation 1: Molecular Dynamics Using Standard AMBEIUAA Charges (First Iteration Charges). Because our main interest in the initial stage was the phosphate-Na+ geometry, we defer a discussion of structural perturbations in the DNA alone until later. Examination of phosphate-Na+ interactions in snapshots from the dynamics trajectory indicated a tendency for the counterions to move away from their initial position along the P-OAR-OB bisector to become more closely associated with just one of the two phosphate oxygen atoms. As the OA-P-OB angle is generally in the range 115.5"-118", it was convenient to produce an initial division of the phosphate-Na+ geometries into those with P-Na+ vectors within 30" of the bisector or, if not, within 30" of either the P-OA or P-OB bond vectors. More diffuse arrangements were categorized having P-Na+ vectors within 50" of either of the phosphorus-ester oxygen bond vectors. Examination of distribution plots (Figure 3) revealed two maxima within 6 8, separated by a minimum at 4 A. These P-Na+ distances were therefore chosen to further subdivide the interaction geometry categories (Table 1). The changes in the populations of each family of conformations over the MD simulation are shown in Table 2. As placed initially by AMBER,35the phosphatekounterion geometries all correspond to family 1A. However, as the simulation progresses the population of this family falls rapidly. There is a clear tendency for the counterions to associate with just one or other of the phosphate oxygens, rather than sit between them. The distribution is not symmetrical, association
Laughton et al.
11594 J. Phys. Chem., Vol. 99, No. 29, 1995 Y."
1 Family I
6.0
sM
1
M
4.0{
2,0u '
0.0
A h
30.0
OD
20.0
I \
4
10.0
30.0 25.0
A
Family 3
,
,
10
20
5.0
0.0
.-.", h
75.0
O9
50.0
25.0 0.0
\
o
as before (see Table 4), the dynamics was again restarted from the point at which the counterion restraints were first removed. After 80 ps (30 ps of partially restrained dynamics, followed by 50 ps of free MD), the populations of the different counterion/ phosphate geometry families were very similar to those observed in the previous simulation, so the MD was continued for a further 100 ps (simulation 3A). In order to obtain more data on the phosphate-Na+ system, a divergent MD simulation was performed (simulation 3B). For this simulation the starting coordinates of the system were taken from the 30 ps time point of the previous simulation (the point at which all restraints were removed) and so were the velocities, except that all the signs were inverted. The change in the sign of the velocities introduces a large perturbation in the system, which would lead to a dramatically different trajectory, but does not introduce changes in the kinetic energy of the system, which avoids problems of thermal equilibration. This second trajectory (referred to as the divergent trajectory) was also followed for 150 ps. The family populations, averaged over corresponding time periods of the two simulations, are shown in Table 5 As can be seen from Table 5, there is little evidence for further drift in the average populations of the different families ( c o n f i i n g the reasonable convergence of the simulation with respect to the ion distribution), but some show considerable (and slow) fluctuations. The average P-Na+ separation also appears equilibrated after 80 ps. From the distribution plots shown in Figure 4 it can be seen that the average configurations of the phosphate-Na+ system are similar in these two MD simulations. The main exception to this is family 3; from Figure 4 it can be clearly seen that the maximum in the radial distribution function (rdf) at about 3.2 8, is twice as high in simulation 3B as in 3A. As noted before, two maxima appear in the distribution plot: one at a distance P-0 between 3.2 and 3.5 8, while the second, more diffuse (see Figure 4)at a distance around 5 8,. Inspection of structures reveals that the first maximum corresponds to a direct phosphate-Na+ interaction, while the second is a solventseparated contact. The minimum appearing in all the configurations at about 4.0 8, corresponds to an unstable situation, where the layer of water separating phosphate and Na+ in the 5.0 8, maximum has disappeared, but the direct phosphate-Na+ interaction geometry is still not optimum. The Na+ is rarely found along the 0-P-0 bisector and interacts more often with OA than with OB. Distribution plots (see Figure 4) show that the configuration of the phosphateNa+ system with highest probability corresponds to that with a direct Na+-OA contact, Le., the first maximum of family 2. As expected, such a preference changes when looking at population factors, since then entropic considerations make families defined using a larger volume of the configurational space more popular. Averaged population values (see Table 5) indicate that families 7, 6B, and 2A are the most popular ones, followed by families 3A and 4A. Configurations with a very small probability are those corresponding to families 6A, lB, 4B, 3B, and 1A. A more detailed inspection of Table 5 reveals that the Na+ is coordinated directly with the OA atom around 18% of the time, but only 11.4% of the time with OB. The difference is larger if only the more lineal arrangements are considered (15% vs 7.9% for families 2 vs 3). All these differences agree well with the corresponding distribution functions, which have a maximum of around 200 for family 2A and around 75 for family 3A. Translation of these number to free energy differences is straightforward from eqs 4 and 5. Thus, the free energy change due to the movement of a Naf from the optimum position where it is directly coordinated to OB to the optimum position for the
0l 0 0 0
30
40
50
60
r (A)
Figure 3. Radial distribution functions for P-Na+ obtained using standard AMBER charges, subdivided by conformational family.
with OA seems preferred over that with OB (even though at this point the charges on both atoms are equal), and when associated with OA a direct interaction (family 2A) seems more populated than a solvent-mediated one (family 2B). For interaction with OB the reverse is true; families 3B and 5B are more popular than families 3A and 5A. By the end of the simulation period over 20% of the phosphatelcounterion geometries correspond to a diffuse and solvent-separated family (6B),and over 30% of phosphates have no counterion within 6 8, of the phosphorus atom at all (family 7). Simulation 2: Molecular Dynamics Using Second Iteration Modified Charges. After recalculating the partial charges for the phosphate/counterion system, the MD simulation was run again (simulation 2), starting from the point at which the restraints on the counterions were first removed. The analysis of the resulting geometries is shown in Table 3. It can be seen that the tendency of the counterions to diffuse more than 6 8, from a phosphate group was much reduced, as was the percentage of loosely-defined interaction geometries (see for instance the ravein Table 3). As before, close interaction between a counterion and OA is more common than interaction with OB, but now the solvent-separated geometries are less populated than those with direct phosphate-counterion interactions for both OA and OB. Simulations 3A and 3B: Molecular Dynamics Using Third Iteration Modified Charges. After recalculating partial charges
. I . Phys. Chem., Vol. 99, No. 29, 1995 11595
Counterion Distribution around DNA
TABLE 2: Populations (as Percentages) of the Different Families during the MD Simulation with AMBEWAA Charges (Simulation 1); Average P-Na Separations Are Also Displayed family
time (ps)
1A
1B
2A
0-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80
99.9 15.2 4.7 6.4 5.7 0.9 0.3 0.9
0.1 17.5 6.5 5.8 0.9 2.8 4.5 7.2
0.0 14.3 14.7 17.3 17.5 18.6 12.3 9.7
2B 0.0 3.2 3.8 3.6 5.5 6.3 7.9 8.0
3A
3B
4A
0.0 5.5 9.5 10.4 9.2 2.1 1.4 2.8
0.0
0.0 10.4 13.6 10.6 2.9 0.0 0.0 0.0
1.7 2.5 3.0 7.8 6.9 14.1 10.0
4B 0.0 0.6 2.4 1.1 0.8 0.4 0.6 0.9
5A
5B
6A
6B
0.0 5.6 5.5 4.2 1.0
0.0
0.0
0.0
4.5 4.3 4.2 4.7 3.3 5.2 5.2
0.3 0.5 0.3 1.8 0.5 0.1 0.2
17.1 19.7 13.5 22.9 27.9 22.6 21.4
0.0 0.0 0.0
7 0.0 4.1 12.1 19.7 19.2 30.4 31.0 33.7
(4
ravg(l
3.10 4.39 4.71 4.81 4.99 5.45 5.74 5.84
Average distance between sodium ion and closest P atom. TABLE 3: Populations (as Percentages) of the Different Families during the MD Simulation with Second-Iteration Charges (Simulation 2): Average P-Na Separations Are Also Displayed family time (ps) 0-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80 80-90
1A 99.9 25.4 10.1 12.8 4.9 4.3 3.4 3.6 3.9
1B
2A
2B
3A
3B
0.1 21.2 3.4 4.2 5.5 6.2 2.4 6.6 8.9
0.0 17.7 18.5 18.0 26.1 26.4 21.1 15.2 26.7
0.0 0.6 3.5 2.5 2.4 3.8 3.9 5.8 0.3
0.0 10.5 11.0 11.5 15.5 18.1 18.2 14.6 12.0
0.0 1.1 3.7 5.0 4.3 3.0 6.8 5.1 3.2
4A 0.0 10.1 9.8 3.2 2.9 2.3 3.6 5.2 6.3
4B
5A
5B
0.0 2.9 2.1 1.4 4.9 1.5 0.5 0.5 1.2
0.0
0.0
1.9 5.8 6.0 4.7 4.9 4.5 3.3 2.9
2.6 3.3 6.4 2.8 3.7 3.3 4.8 2.6
6A 0.0 0.3 1.1 0.4 0.9 0.6 2.1 1.2
0.0
6B 0.0 4.2 22.3 16.4 16.6 12.0 20.6 13.7 12.0
7
0.0 1.5 5.5 12.4 8.5 13.3 9.7 20.4 20.0
rav/ (A) 3.70 4.05 4.40 4.68 4.5 1 4.59 4.57 5.01 5.08
Average distance between sodium ion and closest P atom. TABLE 4: Standard AMBER and Modified (Third Iteration) Partial Charges for the PhosphatdCounterion Svstem atom standard charge modified charge ~
C5' 05'
P OA OB 03' C3' Na+
0.188 -0.509 1.385 -0.847 -0.847 -0.509 0.233 1.000
0.273 -0.591 1.631 - 1.024 -0.992 -0.540 0.279 0.989
direct coordination to OA (eq 4) is -0.58 kcdmol. Furthermore, the transfer of the Na+ from OB to OA coordination (for a range of distances up to 4 A; eq 5) is associated with a modest change of free energy of -0.39 kcdmol. Another interesting result arises from the large population of family 6 (around 15%). This family contains configurations where the Na+ is closer than 6 to the phosphate, but it is not placed along the usual P-OA, P-OB, or OA-P-OB vectors. The larger popularity of subfamily 6B versus 6A indicates that the Na+ in this family does not commonly involve a direct interaction with the phosphate. To see whether interaction with some other atom was commonly associated with this group, the closest DNA atom to each counterion in each MD snapshot was calculated. The results are shown as a histogram in Figure 5. It can be seen that, as expected, OA is the DNA atom most commonly closest to a counterion, followed by OB. Otherwise, the counterion is most likely to be closest to one of the hydrogen atoms in the same region of space, i.e., those associated with C3', C4', and C5'. Examination of counterion density plots (Figure 6) confirms that there is no evidence for family 6 commonly being associated with the interaction of the counterion with a particular DNA atom; Le., the family is truly diffuse. Population factors can also be used a priori to obtain information on the free energy of association between phosphate and the sodium cation. The practical problems are the same as when this is obtained by integration of potential of mean force
profiles, i.e., the definition of the cutoff dividing "associated" and "dissociated" species. If the minimum in the distribution function is considered as a "cutoff" (4A), it can be concluded that around 40-50% of the Na+ are directly associated with a phosphate. However, if the cutoff is set at 6 A, around 7080% of the Na+ can be considered bound to the DNA. Our MD results can be directly compared with previous theoretical data. Thus, MD simulations by Miaskiewicz et al.I6 found around 50% of the sodium ions in direct contact with the DNA, in good agreement with our results, although Miaskiewicz et al. did not use periodic boundary conditions in their simulations. Results can be also qualitatively compared with theoretical methods that use a macroscopic description of the DNA/water/Na+ system. Thus, Manning's the0ry"~9~~ predicts around 75% of the sodium cations to be bound to the DNA, and Fixman4' using a Poisson-Boltzman method suggested that 50%-66% of the Na+ were directly bound to DNA. General Characteristics of the Molecular Dynamics Simulations. The two 150 ps MD trajectories 3A and 3B provide information not only on the phosphate-Na+ system but also on the structure of the whole DNA system during the dynamics. At this point it should be noted that there are few long (> 100 ps) MD simulations of fully hydrated DNA structures published in the literature. To our knowledge, only two MD simulations have been published that involve a hydrated DNA dodecamer: one by Beveridge's g r o ~ p ' ~using 9 ~ ~the GROMOS(SPC) force field4*and restraints to maintain the DNA structure and the other a 150 ps MD simulation by Weinstein and co-workersI6 using the AMBEWAA(TIP3P) force field22v35on a nonperiodic system. In order to complete the analysis and to check the effect of the change in the phosphate charges on the general structure of the DNA obtained during the dynamics, we extended our original simulation (simulation 1) performed using the standard AMBEWAA charges to 150 ps. A detailed analysis of the differences between the DNA trajectories lies outside the scope of this paper, but some general observations are worth making. In Figure 7a are plotted the root-mean-square (rms) deviations of the MD structures from canonical B-DNA as a function of
Laughton et al.
11596 J. Phys. Chem., Vol. 99, No. 29, I995
TABLE 5: Populations (as Percentages) of the Different Families during the MD Simulations with Third-Iterations Charges, Averaged over the Two Traiectories (3A and 3B); Average P-Na Radii Are Also Displayed family time(ps) 1A 1B 2A 2B 3A 3B 4A 4B 5A 5B 6A 6B 7 rav; (4 ~
30-40 40-50 50-60 60-70 70-80 80-90 90-100 100-1 10 110-120 120- 130 130- 140 140-150 150-160 160-170 170- 180 (I
6.3 3.3 3.3 3.2 5.7 6.9 2.3 1.5 5.3 2.8 2.7 2.9 5.3 4.1 5.3
8.9 10.5 5.4 5.6 1.5 3.2 3.5 3.3 5.9 2.8 2.9 0.8 0.6 4.8 3.5
14.8 17.1 20.3 16.5 14.8 12.9 12.9 12.6 12.5 9.8 14.1 15.8 12.9 16.3 16.3
2.5 5.0 2.1 2.0 2.5 1.7 2.3 3.8 3.5 3.5 1.1 2.3 4.2 2.7 4.3
12.6 10.3 9.3 11.2 9.1 6.8 5.7 6.5 4.6 10.6 8.2 9.1 10.8 5.7 5.6
6.1 5.3 5.3 7.9 8.5 6.3 6.1 5.2 4.5 5.6 8.3 9.4 7.4 6.5 5.5
7.6 8.2 5.9 3.9 4.7 4.0 3.7 5.8 6.1 7.1 5.9 2.9 1.6 2.5 4.7
5.2 3.4 4.2 4.3 4.0 5.6 6.2 5.3 4.5 5.2 5.7 7.5 5.9 5.8 5.0
1.8 2.0 2.5 2.6 2.5 1.9 2.0 1.9 2.3 4.1 2.9 2.5 3.5 3.2 2.5
6.5 3.4 6.3 7.0 5.1 4.9 6.4 7.9 5.5 6.3 5.1 4.9 6.2 8.6 8.1
1.2 1.0 1.8 2.4 2.3 1.0 2.3 1.6 1.5 2.3 2.2 2.0 1.5
16.5 18.5 19.4 18.4 21.0 18.8 23.8 24.5 21.8 17.2 21.4 14.8 15.4 16.6 16.7
0.8 1.7
10.3 11.9 14.1 15.2 18.3 26.3 23.0 20.4 22.0 22.8 19.5 25.3 24.6 22.4 20.8
4.66 4.82 4.84 4.88 5.03 5.20 5.22 5.19 5.20 5.09 4.97 5.02 5.05 5.12 5.10
Average distance between sodium ion and closest P atom.
30-
average Ibl) OD
100.0 50 0
0.0 100.0
OA
CB
05' H5'1 H5'2 H
