J. Phys. Chem. 1995,99, 10373-10382
10373
Monte Carlo Simulation Study of Ion Distribution and Osmotic Pressure in Hexagonally Oriented DNA Alexander P. Lyubartsevt and Lars Nordenskiiild* Division of Physical Chemistry, Arrhenius Laboratory, Stockholm University, S-106 91 Stockholm, Sweden Received: October 14, 1994; In Final Form: January 17, 1995@
The electrostatic osmotic pressure in a system of hexagonally packed DNA molecules has been calculated with the Monte Carlo (MC) simulation method. The DNA molecules were modeled as hard cylinders with charged groups located at the sites corresponding to B-DNA, with the ions considered as point charges with repulsive r-I2 potentials and the solvent treated as a dielectric medium. Periodic boundary conditions for a hexagonal cell were used with Ewald summation of the electrostatic interactions. The pressure was calculated from the relation P = -AF/AV, where differences of free energies F were obtained with the expanded ensemble method. The calculations were carried out both for salt-free solutions and for solutions containing added salt; in the latter case the simulations were performed within the grand canonical ensemble. In the system with only monovalent ions, the forces between DNA were found to be always repulsive. In the case of divalent counterions, an effective attraction between DNA molecules may appear for distances of 5-15 8, between the surfaces, depending on the ion size and salt concentration. The results of the simulations showed that a correct statistical-mechanical treatment of the electrostatic interactions in the frame of a continuum dielectric model can reproduce the essential features of available experimental data, indicating that this contribution to the pressure is an important contributor to the experimentally observed pressure versus distance curves for ordered DNA.
1. Introduction Recently, much interest has been focused on the interactions between ordered DNA molecules in different systems, and a number of investigations'-6 related to this problem have been published. The understanding of the stability of such forms is of considerable practical interest, since in living systems DNA is often arranged in highly ordered and compact forms. Many of the physical properties of DNA in solution are strongly dependent on electrostatic interactions between DNA and its surrounding^.'-^ The role and importance of these polyelectrolyte properties, which generally appear as pronounced salt and counterion valency effects, have been less studied in the case of ordered and compact DNA systems. However, it is well documented that an addition of multivalent ions may lead to condensation of DNA into an ordered This may be taken as an indication of a physical mechanism involving the electrostatic interaction between the counterions and DNA,3,4 although different explanations have been p r o p o ~ e d . ~ , ~ * ~ Usually, the electrostatic interaction between highly charged polyelectrolytes is treated within the framework of the classical electrostatic double-layer theory based on the mean-field Poisson-Boltzmann (PB) equation.I0 This theory always gives a repulsive force between uniformly charged macromolecules and is thus not capable of explaining the condensation to an ordered DNA phase induced by multivalent ions. An attractive contribution in such an approach is usually introduced by adding a van der Waals interaction (this is the case of the DLVO theory of colloidal stability"), but this contribution is not expected to depend on counterion concentration or valency and can therefore not explain the appearance of condensed DNA systems. In a number of contributions, Parsegian and co-workers have studied the DNA-DNA interactions in condensed gels of 'On leave from The Scientific Research Institute of Physics, St. Petersburg State University, 198904, St. Petersburg, Russia. @Abstractpublished in Advance ACS Abstracts, June 1, 1995.
0022-365419512099-10373$09.00/0
ordered DNA by means of an osmotic stress t e ~ h n i q u e , * ~ ~ ~ ~ - ~ ~ measuring the force as a function of distance in the systems. On the basis of the fact that these results cannot be reproduced by traditional double-layer theory, the conclusion of these studies was that electrostatic forces did not play major a role in the experimental behavior. Instead, other contributions to the force between DNA molecules, such as hydration and electrostatically mediated configurational fluctuation forces, have been suggested.2ss.6 However, more accurate treatments of the electrostatic interactions, using both the hypemetted chain (HNC) approximationI3-l6 and Monte Carlo simulation^,^,^.'^,'^ have shown that the inclusion of ion correlations that are neglected in the mean-field PB theory gives an attractive contribution to the electrostatic force which in some cases (e.g., for divalent ions) may lead to a net attractive force. The physical origin of this contribution is clear and has a similar background as the quantum mechanical attractive van der Waals interaction between two hydrogen atoms, caused by the correlated interaction between the electron clouds. In the present case, it is the correlation in the interactions between the ion clouds associated with the charged surfaces that gives rise to the attractive Contribution to the force. Thus, the contribution from this physical mechanism must be taken into account in any discussion of forces in experimental systems. Monte Carlo simulations performed in this laboratory3s4have also shown the importance of including correlation effects in the treatment of the electrostatic interactions for describing forces in ordered DNA systems. In order to provide a sound basis for a discussion of the physical origin of the behavior of ordered DNA systems and the importance of hydration and other force mechanisms, we believe that it is necessary to further investigate the behavior and implications of the electrostatic interactions as obtained from a less approximate treatment than the PB model, taking the effects of ion correlations into account. To achieve this goal it is also
0 1995 American Chemical Society
Lyubartsev and Nordenskiold
10374 J. Phys. Chem., Vol. 99, No. 25, 1995
important to investigate methods of calculating the electrostatic force and also the approximations that are involved within treatments going beyond the mean-field PB approach. The effective force between DNA molecules in an ordered array is directly related to the osmotic pressures3 The standard way of evaluating the pressure in computer simulations is based on the vinal equationt9
0
0
0
0
0
0
0
Figure 1. Hexagonally ordered DNA system with Monte Carlo cells
containing different numbers of DNA. where U is the potential energy of the system and ru denote interparticle distances. The problem of using formula (1) in computer simulations is that fluctuations in the second term usually are large and may in fact be greater than the value of the pressure itself. Another way of evaluating the pressure uses the thermodynamic relation
p = - - aF
av
where F is the free energy of the system. To calculate the pressure, free energy differences of the system at different volumes must be evaluated. In the last few years considerable progress has been made in free energy calculation techniques within computer simulation^.^^-^^ Calculations of the pressure from computer simulations as a derivative of the free energy were, e.g., recently made for hard spheresz5and for a LennardJones system.26 In the present work we have used the expanded ensemble method23 for calculation of the free energy differences. In previous studies this method was successfully applied to free energy calculations for electrolyte solutions23and to a LennardJones system and water.24 The advantages of the expanded ensemble method for pressure calculations will be discussed below. The aim of our work is to use the suggested method for pressure calculations of an ordered DNA system with a strict treatment of electrostatic interactions on the basis of the Ewald summation method, taking the long range of the electrostatic interactions into account, and to apply this to model experimental pressure versus distance curves under different conditions, such as varying salt concentration, counterion type, and counterion valency. In comparison with the previous studies from this l a b ~ r a t o r y ,there ~ , ~ are several differences. We use the expanded ensemble method, which is found to be most suitable for the present problem. We have now taken the long range of the electrostatic interactions into account by means of the Ewald summation method. Furthermore, and most important, the applications of the present calculations are directly relevant to the experimental pressure curves, since we have performed the simulations for added salt and at different salt contents. Consequently, it will be possible to discuss the experimental results in relation to the importance of the electrostatic contribution to the osmotic pressure and to the approximations within our model. In order to take effects of ion type to some extent into account, we have also performed calculations for varying finite ion size. The paper is organized in the following way: In section 2 a detailed description of the model is given. Section 3 describes an application of the expanded ensemble method to the pressure calculation. The results of the simulations and a discussion including comparison with experimental data are given in section 4. The conclusions are stated in section 5.
2. The Model
The system that we want to model consists of an assembly of oriented DNA molecules in hexagonal packing. The model that we consider in the present study consist of a number of DNA molecules packed in parallel (along the Z axis) and separated from each other by a distance R, with mobile ions between them. In the perpendicular plane, the DNA molecules form a hexagonal lattice (see Figure 1). Each DNA is treated as an infinitely long hard cylinder of radius a = 9 A with charged “phosphate groups” situated on the surface of the cylinder at the sites corresponding to the B-form of DNA. Each phosphate group has a charge of -e and a soft repulsive r-I2 potential with effective radius a42 = 2 A. This set of phosphate groups then forms an approximate model of double-stranded DNA with two grooves. Due to the repulsive short-range interaction, the average effective radius of DNA is about 10 A. The ions are treated as point charges with a repulsive r-I2 potential of effective diameter ui. The solvent is considered as a dielectric medium with permitivitty c = 78. All simulations have been carried out at a temperature of 300 K. The total potential energy of the system is
+
where bd = kT((ui uj)/2)’*and rij = Iri - rjl. At this stage it is appropriate to comment on the approximation within this model of treating the whole system, Le., the solvent as well as DNA, with a continuous dielectric constant appropriate to bulk water. There are several levels of description of ionic and polyionic solutions. The first level originates from the electrolyte approximation of Debye and Hiickel and uses the full (nonlinearized) Poisson-Boltzmann mean-field theory.I0 At this level the solvent is treated as a continuum and so is the solute in terms of the ion densities and electrostatic potentials. This level is valid mainly in the limit of weak electrolytes at low c ~ n c e n t r a t i o n . ’ At ~ , ~the ~ next level, one can rigourously describe ion-ion interparticle interactions via a dielectric continuum. The simplest case is the primitive electrolyte and our potential as defined in eq 2 is a small modification of that model. The ion-ion and ion-phosphate potentials of this model are approximations to the exact solventmediated effective potentials between the ions in the solution, which in fact represents the free energy of ion-solvent interactions. From this point of view, c in eq 2 can be regarded as a parameter of this effective potential. This 6 does not represent a local dielectric constant, which, according to some estimation^,^^ may assume much lower values (6-30) in the first 1-2 molecular layers near a charged surface, compared to the bulk value of about 78. The ion radii ai have the sense of
Ion Distribution and Osmotic Pressure in DNA
J. Phys. Chem., Vol. 99, No. 25, 1995 10375
implies two types of transitions, particle displacements and hydrated ion radii and are adjustable parameters of the model changes of an expansion parameter, which gives probability that include effects of the molecular nature of the solvent in an distributions over subensembles that are proportional to their averaged form. The most suitable values of ai for specific ions partition functions. The ratio of the partition functions of the can be obtained, e.g., from fitting to osmotic coefficient data subensembles then yields the free energy difference. (see refs 29 and 31). Finally, the most detailed description is In the present case, in order to evaluate the pressure we need an explicit treatment of the solvent at a molecular level. This to calculate the free energy difference of the same system at level, with an appropriate choice of the atom-atom potentials, different volumes. The partition function in the expanded gives a full description of the system within the frame of a ensable can be written as23 classical (nonquantum) theory and may also give a reliable estimation of effective ion-ion potentials. Existing studies on simulation^^^^^^ or liquid-state integral equation t h e o r i e ~for ~~,~~ ion-water systems at the molecular level have shown that the effective (mean force) potentials of ions usually have one to two oscillations around the potential of the primitive model (with (3) E = 78) at small (few angstroms) distances between the ions. These oscillations reflect the molecular nature of the solvent. where {V,} is a sequential set of volumes, Z(V,) is the Still, for many applications, including strong polyelectrolytes configurational partition functions at volume V,, and { q m } is a and high salt concentrations, continuum dielectric models (the set of balancing factors for equilibration of the probability primitive electrolyte model) have been successful in the distribution over subensembles. The final result of the calculadescription of, e.g., osmotic coefficients for and tion of the free energy does not depend on the specific values divalent3I ions in the whole concentration range and the B-Z of these balancing factors; the procedure of their optimal choice transition in DNA35 (see also some recent r e v i e w ~ ~ ~Our . ~ ~ ) . is discussed in detail in ref 23. choice of a continuum dielectric model and the potential in eq In the course of the MC procedure two types of steps are 2 has the intention to separate purely electrostatic effects from made: particle displacements and changes of volume. As a effects of hydration and the molecular structure of the solvent. result we obtain probability distributions pm as a function of volumes, which gives relative free energies In our simulations the MC cell has the form of an hexagonal prism with side L and height h, and it contains several DNA molecules. We have applied periodic boundary conditions in F(V,) - F(V,,) = all directions. Figure 1 shows how to position one, three, four and seven DNA’s in the hexagonal cell in such way that the periodic boundary conditions generate an infinite triangular Taking the numerical derivative of eq 4 with respect to lattice. In the computer program we actually used a paralvolume, we can obtain the pressure in the whole range of lelepipedical cell, formed by the basis vectors el = (LJ3,0,0), volumes { V,} . e2 = (L2/3/2,3L/2,0) and e3 = (O,O,h), which, because of the Equation 3 can, as in the case of the NFT ensemble, be periodic boundary conditions, is equivalent to an hexagonal cell. rewritten so that the integrals at each “m” are taken over the There are several ways of treating long-range electrostatic same (reduced) volume: forces in computer simulations. The simplest one is the minimum-image convention, but this is not sufficient for polyelectrolyte systems, where an inhomogeneous charge distribution exists. More sophisticated approaches take this into account by incorporating a mean electrostatic potential from the neighboring cells,27 while fluctuations of this charge where the transformation xm = xm(x’) transposes the reduced distribution are usually neglected. Another way is summing volume Y into volume V,. For example, for a uniform up the electrostatic interactions from images of charges in the expansion x,,, = x‘V“‘~ we have neighboring cells: but that method is extremely time consuming. In the present work we have used the Ewald summation methodI9 for treating electrostatic interactions. The Ewald sum provides an explicit calculation of the electrostatic interactions from all the charges in the infinite periodic lattice. The special rearrangement of the reciprocal part of the Ewald sum makes it The practical use of eq 5 or 5a in MC simulations is that in possible to greatly reduce the number of operations, making “volume transitions” the particles are considered to remain fixed the Ewald summation even faster for larger systems (with N 2 in space, but parameters of interactions are changed due to 1000) than the minimum-image approach. The details of the substitution of variables x x’ and the term N In V, is added algorithm of the Ewald summation for an hexagonal cell are to the Hamiltonian. On the other hand, if in the course of the given in Appendix 1. volume transition we really move the particles so that they occupy another volume, we have from the detailed balance principle to add the term (V,/V,,)N to the transition probability, 3. Calculation of the Pressure by the Expanded which leads to the same result. In practical calculations the Ensemble Method term N In V, is added to the balancing factors qtm= q, f N The expanded ensemble method for free energy calculations In V,. This term gives a contribution to the pressure P i d = a(N has been described in detail e l s e ~ h e r e ? ~In. ~general ~ the In V,)laV = NkTN,, corresponding to the pressure of an ideal partition function of the expanded ensemble is the sum of gas. partition functions of canonical ensembles for different paramIn the case of the DNA system described in the previous eters, e.g., temperature, volume, or some other parameter of section, the uniform scaling of ion coordinates is impractical, since the hard core of DNA is not changed. The detailed the Hamiltonian. An MC procedure in the expanded ensemble
-
Lyubartsev and Nordenskiold
10376 J. Phys. Chem., Vol. 99, No. 25, I995 description of volume transitions within the present model is given in Appendix 2. The MC procedure in the volume-expanded ensemble resembles that in the NPT ensemble. If we put q m = -PVmPex in eq 3 and replace by JdV, we obtain the exact expression for the configurational partition function of the NET ensemble. For such a choice of qm, the system will fluctuate around a volume V corresponding to given pressure P,,. In general we can write l;lm = -pVmPex(V), which means that we simulate the system under an external pressure, which has a specific dependence on volume. Thus in the volume-expanded ensemble the balancing factors l;lm have a clear physical meaning: they are directly related to the external pressure. As noted in ref 23, the optimum choice of balancing factors is
This corresponds to uniform distribution of probabilities pm.It also means
or
Le., the internal pressure of the system is equilibrated by the external pressure. The simulation procedure in this case is similar to that of the force balance method.25 In general the expanded ensemble method gives free energies and pressures over a range of volumes. However, for calculation of the pressure at a given volume, only the free energy difference between a pair of nearly equal volumes, VOand V I = VO dV, is needed
+
1 PI -1n-+ P A V Po
P,, (7)
where P,, = (70 - l;ll)/pAV. The expanded ensemble here consists only of two canonical ensembles with volumes differing only by a small element, V and V AV. Transitions from one subensemble to another will in this case occur especially easily, and the precise choice of P,, is not crucial, since for a sufficiently small AV, the transitions between the volumes V and V AV will always take place. However, in order to obtain as high accuracy as possible, P,, should be chosen close to the true internal pressure. This is most conveniently done iteratively. The optimized value of Pex then gives the main contribution to the pressure, while the term In@&,), calculated in the MC procedure, only gives a minor correction to the final value. Since in the present approach we calculate free energy differences between nearly identical systems, it is appropriate to ask why the pressure is then not evaluated by some other approach, as, e.g., the perturbation method.21 The point is that if we apply a perturbation expression to the free energy difference between the volumes V and V AV, we obtain exactly (for AV 0) the virial expression (1) for the pressure, and this may give large statistical fluctuations in the evaluation. In practical pressure calculations the expanded ensemble method can be used either over a wide range of volumes or by calculating the free energy difference only between two nearly
+
+
-
+
equal volumes. The first method is more suitable for a detailed calculation of the dependence of pressure on volume. It requires more accurate choice of the balancing factors qm (or external pressure P,,(V)) and becomes more time-consuming for systems with a large number of particles. On the other hand, it can cover a significant range in a single computer run. The second method is more suitable for pressure evaluation at one point or a limited number of points. The presented algorithm gives the total pressure of the polyelectrolyte model described in section 2. It corresponds to the osmotic pressure of an ordered phase of a salt-free polyelectrolyte system. If additional counterions as well as coions are present due to an equilibrium with a bulk salt solution, their pressure contribution is also included. Thus to calculate osmotic pressure in the ordered phase of DNA one should subtract from the total pressure the contribution of added salt. These data can be obtained separately from a simulation of a bulk electrolyte model, although the salt concentration in the bulk phase may differ from that in the gel phase. One way of evaluating the osmotic pressure when added salt is present is to carry out the simulation of both the ordered gel and bulk phases in the grand canonical ensemble, at a fixed chemical potential. A similar technique was used in references 16 and 18 for calculation of the osmotic pressure and the effective force between charged surfaces (the simulations in these studies were actually carried out in an NVT ensemble, the chemical potential being calculated by using the particle insertion method). The simulation of the bulk phase yields two quantities: ion concentration corresponding to a given chemical potential and the contribution of added salt to the pressure. The latter can be subtracted from the pressure of the gel phase to obtain the osmotic pressure. This method was used in the present work for evaluation of the pressure in the presence of salt. We can substitute in formulas (3)-(7) the canonical partition function Z by the grand canonical partition function Z and then calculate the dependence of the thermodynamical potential S2 = F - C a PaNa = -kT In(=) on volume are the chemical potentials of ions of a-type). Then the pressure can be obtained as
p=--
ElT$
The algorithm remains in general the same as in the NVT ensemble calculations described above, with the exception that MC steps with deletiodinsertion of ions are added.
4. Results and Discussion (a) Ion Distributions. Most previous studies of theoretical modeling of counterion-DNA interactions have been made for the “primitive model”, in which a single DNA in a cylindrical cell is treated as a uniformly charged hard cylinder and the ions are considered as charged hard spheres!~9,14,27,37,38 This cell model is generally believed to be a reasonable model for DNA in isotropic solution at a finite polyion concentration, but it has also been used to model concentrated ordered DNA syst e m ~ . ~A* common ~ ~ . ~ ~ interest in these kinds of calculations has been to evaluate the ion distribution around DNA and to compare the appearance of this distribution obtained from different theoretical methods and for different approximations within a given model. The only study of the ion distribution obtained from an explicit calculation on an ordered DNA system that we are aware of is the Brownian dynamics simulation work by G ~ l d b r a n d . In ~ ~order to study the effects of our model,
Ion Distribution and Osmotic Pressure in DNA 10.
I
-”
J. Phys. Chem., Vol. 99, No. 25, 1995 10377
.
.
.
.
.
.
.
.
,
5
10
15
20
25
30
35
40
45
I
7
I
1
9
10
I
E
gP
0.01
0.001
1
~
0
50
r (4 Figure 2. Ion density distributions of +1 and -1 ions around DNA in 0.1 M electrolyte, R = 125 A. Shown are solution of the PB equation
11
12 r (A)
13
14
15
Figure 4. Density distributions of 4-1 counterions for salt-free solution, R = 30 A. Notation as in Figure 2.
(PB), MC simulation of the primitive model of a single DNA in a cylindrical cell (PM), and MC simulation of the model described in section 2 with one and seven DNA’s in the cell. 51
I
0 0
2
4
6
8
10
12
14
I
r (A) 0
5
10
15
20
25 r (A)
30
35
40
45
50
Figure 3. Electrostatic potential around DNA in 0.1 M 1: 1 electrolyte solution. Notation as in Figure 2.
particularly regarding the explicit presence of neighboring DNA molecules and of the discrete polyion charge distribution, we have made test simulations to compare the results with those obtained from the cylindrical cell model. The first set of simulations were made for the case of DNA with the addition of a 0.1 M 1:1 electrolyte. We carried out the following calculations: (1) solution of the cylindrical Poisson-Bolmann equation (cell of radius 65.52 A); (2) Monte Carlo simulation in a cylindrical cell of radius 65.52 A, height 136 A with 180 counterions and 100 co-ions of radius 2.12 A (we used here the MC-self-consistent-field method described in ref 38; (3) MC simulation in a system of DNA molecules packed in parallel (seven DNA’s in the cell), with DNA being treated as a uniformly charged hard cylinder; (4) MC simulation of the model described in section 2. Calculations (1)-(3) were made for the primitive model of DNA with a radius of 10 A. In cases (3) and (4) the DNA-DNA distance was 125 A, which corresponds to the same ion and polyion concentrations as in case (2). The results for ion radical distributions relative to the DNA axis are given in Figure 2, and the potential distributions are given in Figure 3. It can be seen that for large distances from DNA, the ion distributions obtained by the different methods are very similar to each other. Case 3 is not displayed since it practically coincides with cases 1 and 2. This is an additional confiiation
Figure 5. Ion-ion (lines 1 and 3) and ion-phosphate group (lines 2 and 4) radial distribution functions. Lines 1 and 2 correspond to one DNA in the cell. 3 and 4 to seven DNA’s in the cell.
to be added to previous studies, which have found that the Poisson-Boltzmann approximation is well behaved for polyelectrolytes with monovalent counter ion^.^,'^^^^ However, it may be noted that for smaller distances there are considerable differences in the ion distributions obtained from the models that describe the DNA charge distribution as helical compared to those that treat it as a smeared out cylindrical charge. Thus, the introduction of a more realistic DNA charge distribution results in a lower counterion concentration a few angstroms outside the DNA charges. For a small DNA-DNA axial separation (R = 30 A, corresponding to IO-A separation between the DNA’s surfaces) all the ions will always, due to the spatial restrictions, occur in the close vicinity of one DNA molecule. In this case the details of the local structure become essential, and it is not surprising that the primitive model and the model with localization of charges on DNA give different results for the ion distribution; see Figure 4. On the other hand, the number of DNA‘s in the cell has only minor influence on the distribution and on the ion-ion and ion-phosphate group correlation functions (Figure 5). The results for one and seven DNA molecules do not differ within the limit of statistical error. One of the reasons for this is that the system is strongly periodic, and the Ewald method does reproduce this periodicity. It can also be noted that when the DNA-DNA distance is short (e.g., 30 A), then the ion distribution in the space between
10378 J. Phys. Chem., Vol. 99,No. 25, 1995
Lyubartsev and Nordenskiold
80 70
10
, $-,
20
25
30
35
$0
45
50
DNA-DNA distance (A)
Figure 6. Osmotic pressure in the ordered DNA system with + 1 counterions. Lines: one DNA in the cell and ion diameter u = 0 (1). u = 1 4 (2), and u = 4 8, (3). Points: seven DNA's in the cell and u=4A.
the DNA's is nearly uniform for our present model. Furthermore, there are no strong correlations in the positions of the ions and the phosphate groups: the maximum of the corresponding correlation function is at about 1.55 8, (see Figure 5). This result is due to the description of the DNA charge distribution with discrete charges positioned at the sites of the phosphate groups. (b) Osmotic Pressure in Salt-Free Solution. The osmotic pressure in the system was calculated for a salt-free solution with mono- and divalent counterions in a range of DNA-DNA distances between 25 and 50 A, and with one and seven DNA's in the cell. In addition, several simulations with varying number of DNA's in the cell were performed at the distance 30 8,. Simulations with one DNA in the cell covered the whole range of distances in a single run. In the case of seven DNA's in the cell we ordinarily used two subensembles with slightly differing volumes and carried out a separate simulation for each distance R. The typical length of the simulations was 107-2 x lo7 MC steps. The convergence of the calculations and the statistical error were determined by averaging over pieces of the Markov chain and calculating the variance of the averages. The typical statistical error in the pressure was about 5 Wa. The results for monovalent counterions are presented in Figure 6 for the following three values of the ion diameter: u = 0 (point charges, as in the case of the previous simulations from this laboratory3s4),1 8, (corresponding to the typical size of "naked" ions without hydrated shell), and u = 4 8, (which corresponds to a typical size of hydrated ions). It can generally be noted that irrespective of details of the model, the DNADNA interaction is repulsive over the whole range of distances. There is no difference between the results for one and seven DNA molecules in the cell. The size of the counterion exerts some influence upon the pressure; decreasing the ion size leads to decreasing pressure, the difference being larger for smaller DNA separations. A similar behavior was observed in the case of the interaction between planar charged surfaces (ref 15) studied with the HNC theory. Results for divalent counterions are displayed in Figure 7. In the simulation for ion diameter u = 4 8, we also carried out calculations with one and seven DNA's in the cell. The pressure difference between these two simulations was noticeable although not large. To clarify that the results are independent of the size of the Monte Carlo cell (Le., the number of DNA molecules in the cell) we made a more detailed study for the
-20 20
,
1 :, &.4'
25
30
35
40
(A) Figure 7. Osmotic pressure in the ordered DNA system with +2 counterions. Lines: one DNA in the cell and ion diameter u = 0 (l), u = 1 8, (2), u = 4 8, (3), u = 5 8, (4), u = 6 8, (5). Points: seven DNA's in the cell and u = 4 8,. DNA-DNA distance
TABLE 1: Osmotic Pressure and Electrostatic Energy Results from Simulation of an Ordered DNA System with Divalent Counterions (a= 4 di) at Different Numbers of DNA in the Cella no. of DNA's h, 8, P,,,, lo5 Pa U,IINkT UdNkT 1 3 4 7 9 12
102 136 136 170 170 204
-0.08 -1.11 - 1.22 - 1.35 - 1.42 -1.29
170 170
-1.15 -2.18
-1.688 -1.684 -1.685 -1.685 -1.684 -1.683
0.131 0.130 0.130 0.130 0.130 0.130
primitive model 7,a=98, 7 , a = 108,
statistical error 0.05 0.001 0.0005 DNA separation R = 30 8,;h is the height of the MC cell. axial DNA-DNA distance set to the value R = 30 8,. The results are given in Table 1. It can be seen that, when the number of DNA molecules in the cell is larger than four, there are only negligible differences in the pressure and the total energy obtained from systems of different size, and no significant trend can be found. At this point it is of relevance to ask how it is possible that simulations with a single DNA in the cell can correctly reproduce forces between DNA's in an infinite system and give an attractive pressure. The reason is that the application of the Ewald summation method effectively simulates an infinite system of hexagonally packed DNA. The only approximation is that of introducing periodic constraints on the positions of the ions in different cells. However, as it follows from the results of our simulations, this does not strongly affect the results. It is mainly the ions in the region between two neigbouring DNA's that give the major contribution to the force between the DNA molecules, and the positions of these ions are not constrained by the periodicity imposed by the Ewald summation. In most of our calculations the sites of the charged groups on different DNA's in a cell were correlated, and the actual periodicity of the fixed charges was the same as in the case of a single DNA in the cell. We carried out five calculations in which each of the seven DNA's in the cell was shifted randomly along the Z-axis (for a given value of the distance R = 30 A). The results were in the range from -1.17 x lo5 to -1.72 x lo5 Pa to be compared with -1.35 x lo5 Pa for unshifted
Ion Distribution and Osmotic Pressure in DNA
J. Phys. Chem., Vol. 99, No. 25, 1995 10379
TABLE 2: Results of Simulation of the Bulk Electrolyte in the Volume-Expanded p V T Ensemble with Two Volumes, VI and Vza V2, nm3 P , io5Pa 1 1 (PM)b cid,MIL VI, nm3 1:l Electrolyte a(+l) = a(-1) = 4 8, 1838 212.8 222 0.1000 0.1002 4.73 0.800 0.793 (0.1002) 0.08 1767 229.8 135.7 141.8 0.5 103 0.5102 25.34 0.784 0.780 (0.495) 0.4 220.8 1.0125 1.0117 55.30 0.889 0.881 (1.003) 269.3 280 229.8 0.9 220.8 2:l Electrolyte a(+2) = a(-1) = 4 8, 167 172.2 0.00981 0.00982 0.66 0.723 0.710 (0.0101) 0.0071 9422 9707 0.0490 0.0490 3.12 0.561 0.559 (0.0477) 293.3 301.2 0.0275 3313 3401 0.428 0.416 (0.259) 392.2 406.6 0.2458 0.2449 15.48 0.105 883.3 919 2: 1 Electrolyte a(+2) = 6 A, a(- 1) = 4 A 235. 24 1 0.0049 1 0.00491 0.37 0.795 0.0039 26500 27212 158. 165. 0.099 0.0992 6.53 0.555 919 0.055 883.3 Cid is the density of an ideal gas corresponding to the same chemical potential /?p = In (cidA3),A = h(2~rmk7")-'/~ (in the program, the chemical potential was introduced through Cid). The indices 1 and 2 refer to the corresponding volume. 1 = exp@(u - pid)) = C , d k is the activity coefficient. The last column shows the activity coefficients for the primitive model of the electrolyte (charged hard spheres of diameter 4.25 8,) from ref 27. * Ion concentration is shown in parentheses.
DNA's. The deviations are larger than the statistical error; nevertheless they are small enough compared to the total changes in the pressure in the range of the interaxial DNADNA distances studied and also compared with the "ideal gas" contribution Pid = NkTlV, which is equal 23.2 x 105 Pa for R = 30 A. For the same system (seven DNA's in the cell and R = 30 A), we also carried out simulations with a smeared out charge distribution model of DNA, these of the osmotic pressure results which are shown in Table 1. The deviations are not large, and the conclusion is that the details of charge distribution on the surface of the polyion do not strongly affect the pressure, although it is important for the counterion distribution around DNA at short DNA-DNA separations. Since the results are practically independent of the number of DNAs in the cell, other simulations with different values of the ion diameter u were carried out with a single DNA. It can be seen that for divalent ions the ion size has a crucial influence on the pressure for distances less than 30 8, (note that for larger DNA separations the pressure in all the considered cases is close to zero). For small ion diameters we obtained a strong attractive force at small DNA separation (such a result was also observed in ref 4). For u = 4 8, we have a region of attraction when the distance between DNA is in the range between 27 and 33 A. Note that this region of negative pressure was obtained in a simulation with seven DNA's in the cell; simulation with a single DNA yields a region of attraction between 28 and 3 1 A. Still, the absolute difference in the pressure is quite small between the systems with one and seven DNA's in the cell and is mainly due to a small value of the pressure in this region. The effect of ion size is much larger, so that increasing the ion diameter from 4 to 5 A (and of course also to 6 A) leads to a disappearance of the attractive region. The experimental results for divalent ions exhibit different behavior of the osmotic pressure depending on the ion type present in the ordered DNA system. Thus, for MgC12 the interaction is always repulsive2 and resembles the results that we obtained for u = 5 or 6 A, given in Figure 7. On the other hand, some divalent ions, e.g., MnZ+have been found to cause a spontaneous DNA assembly into an ordered phase.5 Our result for an ion diameter 0 = 4 8, strongly resembles the experimental osmotic pressure in 0.05 M MnClz solution at 35 "C, including the region of attraction between 27 and 33 8, (of course, there is a gap in the experimental pressure graph in the region of attraction). Available data on the evaluation of the effective ion size from osmotic coefficient^^^ show that the effective Mg2' hydrated radius is indeed larger than that of hydrated Mn2+. A more detailed quantitative comparison with experiments is given in the next section for the case af added salt.
The qualitatively different behavior observed experimentally for the forces between oriented DNA molecules in systems with varying divalent ions can be explained in this approach by size differences in the hydrated shell of the ions. A small change of the ion size just in the range of the typical values of the hydrated ion diameters may cause a transition from a purely repulsive interaction to an attractive one with a spontaneous assembly of DNA into an ordered phase. (c) Added Salt. To study the effect of added salt due to the ordered DNA system being in equilibrium with a bulk electrolyte solution, some simulations were carried out in the volumeexpanded pVT ensemble with addition of 1:1 and 2:1 salts. This kind of simulation in the grand canoncal ensemble is a natural method to use for mimicking the experimental situation where the salt concentration of the ordered DNA phase is not known. Instead it is given by the thermodynamic constraint that the chemical potential in the gel phase must be the same as in the bulk electrolyte phase with which it is in equilibrium and whose concentration is experimentally given. In the simulations, transitions with deletion or insertion of ions were peformed by using the standard procedure (see, e.g., ref 28), adding or removing two or three ions simultaneously to keep the system electroneutral. For all cases we made separate simulations of the ordered DNA system (with one DNA in the cell of height h = 102 A) and of the bulk phase (without DNA) with equal chemical potentials. The effective hydrated ionic diameter u was set equal to 4 8, for both mono- and divalent counterions in the DNA system as well as in the bulk phase. The simulation of the bulk phase at a given chemical potential provides the corresponding ion concentration and gives the contribution of the added salt to the pressure. The results of the simulations of the bulk electrolyte phase for different chemical potentials are given in Table 2. It can be noted that for our ion model with an r-I2 repulsive potential, we obtained activity coefficients very close to those calculated in ref 28 for a primitive electrolyte model. In the results presented below for the osmotic pressure, the contributions of added salt have been subtracted. The results of the osmotic pressure for ordered DNA in equilibrium with an 1:l electrolyte bulk phase are shown in Figure 8 with a logarithmic scale for the pressure. For a nonzero salt concentration the dependencies are nearly linear, which means an exponential decay of the pressure as a function of distance. Our results are in a good agreement with the experimental work of Rau et al.,2 which is also displayed in Figure 8 and in which exponential decay of the osmotic pressure with distance was likewise observed. For an effective force of
Lyubartsev and Nordenskiold
10380 J. Phys. Chem., Vol. 99, No. 25, 1995 le+07 35
-G 3
6
30
Ln
*
0
le+06
25
w
a
W
vl
W
c F 2 100000
2
20
c
15
v)
10
B
v 1 5 0
0 a
-5
YY""
20
25
35 PO DNA-DNA distance (A)
30
45
50
Figure 8. Osmotic pressure in the ordered DNA system in equilibrium with a 1: 1 electrolyte bulk phase at different salt concentrations. Ion diameter is 4 8,. Filled squares are experimental results2 for 0.1 and 0.5 M equilibrating bulk salt solutions.
20
30 DNA-DNA distance (A)
25
40
35
Figure 10. Osmotic pressure in the ordered DNA system in equilibrium with a 2:l electrolyte bu!k phase at different salt concentrations. Diameter of all ions is 4 A. Filled squares are experimental results for 0.05 M M I I C I ~ . ~ 1e 4 7 Q
4
I
4, W
1
le+06
k
25
30 35 -40 DNA-DNA distance (A)
45
50
Figure 9. Average number of co-ions in the MC cell of height h = 102 8,containing one DNA and 1: 1 or 2: 1 electrolyte, as a function of DNA-DNA separation and given for different salt concentrations of the equilibrating bulk electrolyte phase. Ion diameter is 4 8, in all
10000
'
20
I 25
30 DNA-DNA distance
(A)
35
40
cases.
Figure 11. Osmotic pressure in the ordered DNA system in 2:l electrolyte. Ion diameters: o(+l) = 4 8,, a(+2) = 6 8,. Salt concentrations: (0) 0.005 M; (0)0.1 M; filled squares are corresponding experimental results for MgC12 from ref 2. Dashed line is the PB result taken from ref 3.
interaction between DNAf = P,s,R/31'2 in a 0.5 M electrolyte, we obtained an exponential decay constant of 3.4 A, which is the same result as that reported in ref 2. As in the experiment, we obtain larger repulsive forces for low salt concentration. Some differences from the experiments are observed at high concentration (1 M), for which the simulations give a quick fall of the osmotic pressure as a function of distance, with a decay constant of about 2.5 A, whereas the experimental results for 1 M are the same as for 0.5 M. Figure 9 shows the average number of co-ions in the MC cell for different ion concentrations. The cell includes three turns of a DNA helix, so the number of neutralizing counterions is 60 for monovalent and 30 for divalent counterions. One can see that the number of co-ions strongly correlates with the decrease of pressure due to the addition of salt. For short distances between the DNA molecules, the number of co-ions is essentially less than the number of counterions, that is why distribution of ion clouds around DNA and, consequently, the osmotic pressure has only a weak dependence on salt concentration. This dependence becomes stronger at larger separations (of course, for very large separations the osmotic pressure tends to zero in any case).
Results of the 2:l electrolyte simulation are presented in Figures 10 and 11. The hydrated diameter of the monovalent ions was set to 4 A, and for the divalent ions it was 4 A (Figure 10) and 6 8,(Figure 11). Experimental data for DNA in MnC12 (ref 5) and MgC12 (ref 2) are also displayed. In the case of the 2:l electrolyte, the osmotic pressure was less dependent on salt concentration than for the 1:1 electrolyte. The results for 0.005 and 0.01 M practically coincide with the corresponding results (with the same counterion size) for the salt-free solution. This is not surprinsing, since even at the maximum distance of 50 A, the fraction of co-ions is only about 3%. For larger salt concentrations and counterion diameter (5 = 4 A, we have a stronger attractive contribution to the force, especially for inter-DNA distances in the 26-32-A range. It can be noted that practically the same region of attraction was observed in the experiments for ordered DNA in the presence of Mnc11.~ For distances in the range 35-40 A, the results reveal a very weak repulsion, the osmotic pressure being about 30 kPa at 0.05 M and 10-20 kPa at 0.25 M, which, however, is close to the statistical error (about 10 kPa). Results of simulation of the system with divalent ions of diameter o ( 2 f ) = 6 8, are shown in Figure 11 and compared
Ion Distribution and Osmotic Pressure in DNA with experimental results from ref 2 for DNA in MgC12 solution and with results of the Poisson-Boltzmann (mean field) approximation from ref 3. Once again one can see that the MC results reproduce the experimental situation much better than the PB theory.
5. Conclusions This study has been able to reproduce qualitatively, and to some extent also quantitatively, the main typical features of the interaction of charged DNA macromolecules summarized in ref 5: an exponential decay of the osmotic pressure with distance including similarity of the decay constant; the weak dependence of the .results on salt concentration at small distances and stronger dependence at larger distances; a possibility (depending on ion size) of attraction of charged DNA polyelectrolytes in a system with divalent ions. The mean-field treatment of the electrical double layer cannot give an effective attraction (see, e.g., a recent papefl2), but this effect is contained in the present correct statistical mechanical treatment of the electrostatic interactions within the frame of a continuum dielectric model. This fact is a strong argument for the viewpoint that electrostatic forces give an important contribution to this feature of the experimental data. The c o n c l u ~ i o n ~that - ~ the interactions between charged DNA molecules cannot be explained by electrostatic forces, and the consequent need for a contribution of the “hydration force” to explain the main features of the data, is not supported by the present study. However, this does not mean that we exclude the possibility of hydration effects due to the molecular nature of the solvent water. In the previous MC studies from this l a b ~ r a t o r ythe ~ , ~counterions were treated as point charges (a = 0). The present study has shown that this treatment leads to an overestimate of the attractive interaction force for divalent counterions. The inclusion of a hydrated ion radius is a simple way of including the molecular nature of water within the continuum model, which can be regarded as a parameter within our model to account for hydration effects. Further refinements of the model may include a more adequate choice of ion-ion and ion-DNA potentials on the basis of effective (mean force) potentials of ions in water33,34and, as a next step, an explicit treatment of the water molecules. The latter improvement would require much more computer time, although, as follows from recent studies (see, e.g., ref 43 and references therein), simulations of a DNA environment with water and counterions with distances up to 10 A from the DNA surface are possible on modem computers. This paper has also shown a new prospective approach to the pressure calculations based on the expanded ensemble method. The method can be directly applied to an arbitrary statistical-mechanical system, including one with hard-core potentials for which the use of eq 1 is impractical due to the large statistical uncertainty introduced by the difference between two large terms. Actually, a similar approach, only in other terms (using a force balance method), was recently applied to the study of a hard-sphere The present approach uses a more general method for the free energy calculation. Application of the method to the pressure calculation in molecular dynamics simulations is also straightforward in the manner of ref 24, with introduction of a volume scaling instead of temperature scaling.
Acknowledgment. We are grateful to Prof. P. N. VorontsovVelyaminov for helpful discussions. This work has been supported by the Swedish Natural Science Foundation (NFR) and by a postdoctoral fellowship (to A.L.) from the WennerGren Foundations.
J. Phys. Chem., Vol. 99, No. 25, 1995 10381 Appendix 1. Treatment of Electrostatic Interaction in a Hexagonal Cell The periodic lattice of a hexagonal cells is generated by the basis vectors el = (L2/3,0,0), e2 = (L2/3/2,3L/2,0), and e3 = (O,O,h), where L and h are the sizes of an MC cell in the X Y and Z directions, respectively. The Ewald sum can be written as
where a is the Ewald parameter, which is chosen so that the first sum can be taken only over interactions inside rCut= min(L2/3/2,h/2), and the second sum is taken over the reciprocal lattice vector k = nlkl n2kz n3k3, k2 < k2cut,where 1 ~ 1 , 1 1 2 , and n3 are integers (excluding the case nl = n2 = n3 = 0) and the basis vectors in the reciprocal space are
+
2n
+
2x
2n
k, = -e2 x e3, k, = -ije3x e,, k, = -e, V V
x e2
We use a = 3/rcutand kUt = 6a,which gives about 1000 vectors in the k-space K, depending on the ratio uh. The second term of the Ewald sum can be rewritten as
(-42)
This reduces the number of operations from %KiV to %KN for calculating the energy and from xKN to x K for calculating the energy difference after displacement of an ion (in the last case it is necessary to store in the computer memory the values N of C,=,Z, cos(kr,) and ~ ~ ,sin(kr,) z , for each value of k and update them after each successful move; the total number of operations will still be of the order of K ) . If the number of particles N becomes larger than K, the first part of the Ewald sum requires more computer time than the second, in this case one can decrease rcut. Although this leads to some increase of computer time for calculation of the second part, the total time may become less. Appendix 2. Volume Transitions in the System of DNA Molecules Packed in Parallel Consider a transition in which the DNA-DNA distance is changed from R to R’. In principle one can use a uniform expansionkompression of the system in the X Y plane. However, the probability of such transitions would be very low due to the hard core of DNA (the DNA diameter does not change in the volume transitions). Therefore we used a more sophisticated procedure. The whole volume is divided into hexagonal cells around each DNA. Each DNA together with its subcell is moved in a uniform way, so that DNA-DNA distances become equal to R‘. Then all the subcells increase or reduce their borders so that the cells are brought into contact again. The second step is made according to the following rule. Ions inside the cylinder of radius X = R12 change their distance r to the
10382 J. Phys. Chem., Vol. 99, No. 25, 1995
Lyubartsev and Nordenskiold
closest DNA along the line connecting the ion and DNA axis by
/=
[
+
]
(X2-a2)r2 (x2-xt2)a2 (X2-a2)
(A31
(this transformation provides {(a) = a and { ( X ) = X’). The ions outside the cylinder of radius X (at the comers of the hexagon) expand or compress uniformly:
X’
r’ = x r
(A41
The Jacobian of the transformation (A3) is (X‘ - a2)/(XZa2)(note that this is a single transformation where the Jacobian does not depend on r) and the Jacobian of (A4) is (X’/W2.Thus according to ( 5 ) the multiplier ((R’2-a2)/(R2-a2))N1(R’2/R2)N2 should be included with the transition probabilities (N1 is the number of ions inside the cylinders of radius X around each DNA, and N2 is the number of ions outside them).
References and Notes (1) Wilson, R. W.; Bloomfield, V. A. Biochemistry 1979, 18, 2192. (2) Rau, D. C.; Lee, B.; Pagsegian, V. A. Proc. Nari Acad. Sci. USA. 1984, 81, 2621. (3) Guldbrand, L.; Nilsson, L. G.; Nordenskiold, L. J . Chem. Phys. 1986, 85, 6686. (4) Nilsson, L. G.; Guldbrand, L.; Nordenskiold, L. Mol. Phys. 1991, 72, 177. (5) Leikin, S.; Parsegian, V. A,; Rau, D. C.; Rand, R. P. Annu. Rev. Phys. Chem. 1993, 44, 369. (6) Podgomik, R.; Rau, D. C.; Parsegian, V. A. Biophys. J . 1994, 66,
962. (7) Record, M. T.; Mazur, S. J.; Melancon, P.; Roe, J.-H.; Shaner, S. L.; Unger, L. Annu. Rev. Biochem. 1981, 50, 997. (8) Manning, G. S. Quart. Rev. Biophys. 1978, 11, 179. (9) Frank-Kamenetskii, M. D.; Anshelevich, V. V.; Lukashin, A. V. Uspekhi Fiz. Nauk (Moscow) 1987, 151, 595. (10) Reiner, E. S.; Radke, C. J. Adv. Colloid Interjace Sci. 1993, 47, 59. (11) Venvey, E. J. W.; Overbeek, J. T. G. Theory of the stability of lyophobic colloids; Elsevier: Amsterdam, 1948. (12) Podgomik, R.; Rau, D. C.; Parsegian, V. A. Macromolecules 1989, 22, 1780.
(13) Kjellander, R.; Marcelja, S. Chem. Phys. Lett. 1984, 49, 112. (14) Steven, C. L.; Glenn, T. M. Adv. Chem. Phys. 1984, 56, 141. (15) Feller, S. E.; McoQuarrie,D. A. Mol. Phys. 1993, 80, 721. (16) Kjellander, R.; Akesson, T.; Jonsson, B.; Marcelja, S. J . Chem. Phys. 1992, 97, 1424. (17) Guldbrand, L.; Jonsson, B.; Wennerstrom, H.; Linse, P. J . Chem. Phys. 1984, 80, 2221. (18) Bolhuis, P. G.; Akesson, T.; Jonsson, B. J . Chem. Phys. 1993, 98, 8096. (19) Allen, M. P.; Tildesley, D. J. The computer simulation ofliquids, 2nd ed.; Clarendon: Oxford, 1987. (20) Beveridge, D. L.; DiCapua, F. M. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431. (21) Pearlman, D. A.; Kollman, P. A. J. Chem. Phys. 1989, 91, 7831. (22) Strastmaa, T. P.; Berendsen, H. I. C.; Postma, J. P. M. J. Chem. Phys. 1986, 85, 6720. (23) Lyubartsev, A. P.; Martsinovskii, A. A.; Shevkunov, S. V.; Vorontsov-Velyaminov, P. N. J . Chem. Phys. 1992, 96, 1776. (24) Lyubartsev, A. P.; Laaksonen, A.; Vorontsov-Velyaminov, P. N. Mol. Phys. 1994, 82, 455. (25) Attard, P.; Moule, G. A. Mol. Phys. 1993, 78, 943. (26) Valleau, J. P. J. Chem. Phys. 1993, 99, 4718. (27) Murthy, C. S.; Bacquet, R. J.; Rossky, P. J. J . Phys. Chem. 1985, 89, 701. (28) Valleau, J. P.; Cohen, L. K. J . Chem. Phys. 1980, 72, 5935. (29) Triolo, R.; Grigera, R.; Blum, L. J . Phys. Chem. 1976, 80, 1858. (30) Hansen, W. N. J . Electroanal. Chem. 1983, 150, 133. (31) Corti, H. R.; Femandez-Prini, R. J . Chem. Soc., Faraday Trans. 2 1986, 82, 921. (32) Pratt, L. R.; Hummer, G.; Garcia, A. E. Biophys. Chem. 1994,51, 147. (33) Guardia, E.; Rey, R.; Padro, J. A. J . Chem. Phys. 1991, 95, 2823. (34) Pettitt, B. M.; Rossky, P. J. J . Chem. Phys. 1986, 84, 5836. (35) Soumpasis, D. M.; Robert-Nicoid, M.; jovin, T. M. FEBS Lett. 1987, 213, 341. (36) Rashin, A. A.; Bukatin, M. A. Biophys. Chem. 1994, 51, 167. (37) Vlachy, V.; Haymet, A. D. J. J . Chem. Phys. 1986, 84, 5874. (38) Vorontsov-Velyaminov, P. N.; Lyubartsev, A. P. Mol. Simulation 1992, 9, 285. (39) Wennerstrom, H.; Jonsson, B.; Linse, P. J . Chem. Phys. 1982, 76, 4665. (40) Guldbrand, L.; Nordenskiold, L. J . Phys. Chem. 1987, 91, 5714. (41) Guldbrand, L. Mol. Phys. 1989, 67, 217. (42) Overbeek, J. T. G. Mol. Phys. 1993, 80, 685. (43) Swaminathan, S.; Ravishanker, G.; Beveridge, D. L. J . Am. Chem. SOC. 1991, 113, 5027. Jp9427771
