Dielectric Saturation of the Ion Hydration Shell and Interaction

Apr 18, 2007 - The Journal of Physical Chemistry B 0 (proofing), ... Combining molecular dynamics and lattice Boltzmann simulations: a hierarchical co...
0 downloads 0 Views 243KB Size
Download PDF
5264

J. Phys. Chem. B 2007, 111, 5264-5276

Dielectric Saturation of the Ion Hydration Shell and Interaction between Two Double Helices of DNA in Mono- and Multivalent Electrolyte Solutions: Foundations of the E-Modified Poisson-Boltzmann Theory Sergei Gavryushov* Engelhardt Institute of Molecular Biology, 32 VaViloVa Street, Moscow, Russia ReceiVed: October 30, 2006; In Final Form: February 22, 2007

Potentials of mean force between single Na+, Ca2+, and Mg2+ cations and a highly charged spherical macroion in SPC/E water have been determined using molecular dynamics simulations. Results are compared to the electrostatic energy calculations for the primitive polarization model (PPM) of hydrated cations describing the ion hydration shell as a dielectric sphere of low permittivity (Gavryushov, S.; Linse, P. J. Phys. Chem. B 2003, 107, 7135). Parameters of the ion dielectric sphere and radius of the macroion/water dielectric boundary were extracted by means of this comparison to approximate the short-range repulsion of ions near the interface. To explore the counterion distributions around a simplified model of DNA, the obtained PPM parameters for Na+ and Ca2+ have been substituted into the modified Poisson-Boltzmann (MPB) equations derived for the PPM and named the -MPB (-MPB) theory. -MPB results for DNA suggest that such polarization effects are important in the case of 2:1 electrolyte and highly charged macromolecules. The three-dimensional implementation of the -MPB theory was also applied to calculation of the energies of interaction between two parallel macromolecules of DNA in solutions of NaCl and CaCl2. Being compared to results of MPB calculations without the ion polarization effects, it suggests that the ion hydration shell polarization and inhomogeneous solvent permittivity might be essential factors in the experimentally known hydration forces acting between charged macromolecules and bilayers at separations of less than 20 Å between their surfaces.

Introduction Electrostatics of highly charged macromolecules has been the subject of growing interest through the last two decades. Despite reported molecular dynamics (MD) simulations of DNA with water and moveable ions,1-4 all-atom model simulations are still extremely demanding of computer time for many macromolecular systems of interest. Therefore, the implicit solvent description is used in MD5 or Monte Carlo (MC) simulations6,7 to describe ionic solution in contact with charged macromolecules. This level of simplification is also used in all practically applicable theoretical approaches of ion-macromolecule interaction such as the Poisson-Boltzmann (PB) approximation,8,9 modified Poisson-Boltzmann (MPB) theory,10,11 and hypernetted chain (HNC) approximation.12,13 The implicit solvent model treats the solvent as a dielectric continuum. In MC simulations, the dielectric constant has to be homogeneous throughout the system for computational reasons,14 whereas the theoretical approaches allow inhomogeneity of permittivity.15 For example, the interior of a protein molecule is treated as low polarizable dielectric, which is necessary for the PB predictions of pK values of charged groups.16-18 Other dielectric effects such as ion-image interactions are neglected in the PB approximation, but they are naturally described by the MPB equations.11 These interactions are electrostatic repulsion of an ion from the interface between water and nonpolarizable medium. Among the polarization effects on the interactions of macromolecules, there is one that has remained questionable through * Author to whom correspondence should be addressed. E-mail: [email protected].

decades. This is the dielectric saturation of the solvent in the fields of macromolecules and ions. In particular, it leads to the polarization deficiency of the ion hydration shell in the external field, which causes a linear decrease of electrolyte solution permittivity as a function of the electrolyte concentration.19-21 In the vicinity of a highly charged macromolecule such as DNA, the dense cloud of counterions should notably decrease the mean permittivity near DNA and affect the free energy of the system. Experimental dielectric properties of electrolyte solutions have been successfully reproduced by MC simulations of the hydration shell polarization deficiency of Na+, Ca2+, and Cl- in SPC/E water.22 A notable impact of the ion hydration shell polarization on the thermodynamics of an electrolyte solution has been discussed since works by Levine and co-workers.23,24 It was verified by recent MD and MC simulations of interionic potentials of mean force (PMFs) and ion activity coefficients of solutions of MgCl2, CaCl2, SrCl2, and BaCl2.25 On the basis of MC simulations, the primitive polarization model (PPM) of hydrated cations was proposed in ref 22. In addition to the hardcore diameter for ion-ion interactions used in the primitive model (PM) of the electrolyte, this approximation treats ions with their hydration shells as low polarizable dielectric spheres of ion-specific radii (∼4 Å) and of internal permittivities (∼10). This approach was initially called the “dielectric sphere model”.23 Further improvement of the PPM was discussed in ref 26 according to results of MD simulations of interionic PMFs.25 The solvent dielectric saturation effects were approximately incorporated into PB calculations of DNA electrostatics.27 As was mentioned above, the PB approximation is a poor treatment of forces acting on an ion in an inhomogeneous dielectric. An

10.1021/jp067120z CCC: $37.00 © 2007 American Chemical Society Published on Web 04/18/2007

Interaction between Two Double Helices of DNA

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5265

attempt to incorporate an approximate description of the solvent dielectric saturation and polarization of the ion hydration shells into the MPB equations showed that resulting effects on DNA electrostatics could be significant in the case of divalent counterions.22,28 In the present work these polarization effects are systematically explored to be correctly substituted into the MPB equations. The research includes MD simulations of the PMFs between a large spherical macroion and single Na+, Ca2+, and Mg2+ ions in SPC/E water to verify parameters of the PPM and approximate the dielectric model of the ion exclusion zone near the macroion’s surface. It also includes MD simulations of the mean electrostatic force acting on a water molecule around the macroion to verify the implicit solvent dielectric model. Although the first attempts to incorporate the ion polarization effects into the MPB theory were done a half a century ago,23,29-32 the theoretical basis of such calculations with spatial inhomogeneity of the solution permittivity (r) has remained unclear. Therefore, the theoretical foundations of the “-MPB” theory based on the PPM are also considered. Finally, the resulting -MPB equations and PPM parameters are applied to calculate energies of interaction between two double helices of DNA at different interaxial distances and ionic concentrations. Results are compared with those obtained without effects of the ion hydration shell polarization and with the PB calculations. The PPM is still a crude approximation of ions, so the main objective of this study is estimation of a possible impact of the dielectric saturation of the cation hydration shells on the free energy of the macromolecule. It justifies an application of a simplified model of DNA in the present study of DNA-DNA interactions. A comparison of ordinary MPB calculations for PM ions and -MPB calculations for PPM ions could clarify the role of such polarization effects in the ion specificity of electrostatic interactions of highly charged macromolecules. Dielectric Approximation of Ion-Macromolecule Interactions To use the dielectric model of the macromolecule/water interface in PB or MPB calculations, this model should be verified by all-atom MD simulations of water polarization around the macromolecule. The dielectric models of solvent and ion hydration shells can also be verified by a comparison of MD-simulated ion-macromolecule PMFs with electrostatic energies calculated for the implicit solvent model of ionmacromolecule interactions. Such comparisons should be made for different intensities of the electric field of the macromolecule. This analysis is done in the present work for SPC/E water, a model of an 8 Å radius spherical macroion, and Na+, Ca2+, and Mg2+ ions represented by the Lennard-Jones (LJ) approximations of ion-water interactions. In the first part of this section, the MD-simulated PMF acting between a single ion and a charged spherical macroion is compared with the energy from continuum electrostatic calculations. It will allow us to verify parameters of the PPM of the ion and approximate the dielectric discontinuity between the nonpolarizable macroion and water. MD simulations of the water electrostatic potential around the macroion are described in the second part.

Figure 1. Illustration of the model of MD simulations of the mean force between the macroion and the ion at infinite dilution. The system is composed of a spherical cavity containing 2012 water molecules, an ion and a macroion at fixed positions, surrounded by a dielectric medium. The parameter am ) 8 Å means the position of the center of the LJ repulsion from the macroion. R ) 25.0 Å, RRF ) 27.8 Å, dm ) 3 Å, and RF ) SPC/E ) 71. See text for further details.

Figure 2. Illustration of the dielectric model. The value am means the same as in Figure 1, Qm and qi the charges of the macroion and ion, respectively, am the approximating dielectric boundary between water and the interior of the macromolecule, ai and i the radius and dielectric constant of the low permittivity sphere of the hydrated ion, respectively. See text for further details.

Verification of the Continuum Dielectric Approximation for Ion-Macroion Interactions. Models. All-atom MD simulations of ion-macroion PMFs (Figure 1) are compared to calculations of the electrostatic energy of an implicit solvent system (Figure 2). In atomistic model simulations, a macroion, an ion, and 2012 water molecules are enclosed in a spherical cavity of radius RRF surrounded by a dielectric medium with the permittivity RF representing water. Details of this method of PMF evaluation are given in ref 25. In addition to the Coulomb interaction, the macroion interacts with all particles via the short-range water oxygen LJ potential centered at the macroion radius am ) 8 Å. Water molecules interact with the surroundings via the shortrange water-water potential centered at radius RRF ) 27.8 Å. Inside the spherical cavity, all interactions were evaluated explicitly. The LJ parameters (σ and ) of the rigid SPC/E water model33 and ions34-35 are given in Table 1. The LorentzBerthelot mixing rules are applied. The leading reaction field of the dielectric medium acting on water molecules, an ion, and a macroion was represented by image charges (Figure 1)36 using the permittivity of the SPC/E water model, SPC/E ) 71.22,25,37,38 The macroion center was fixed at a distance of dm ) 3 Å from the center of the spherical cavity, whereas the center of

5266 J. Phys. Chem. B, Vol. 111, No. 19, 2007

Gavryushov

TABLE 1: Lennard-Jones Parameters σ and E for Ions and SPC/E Watera atom/ion

σ



O(SPC/E) Na+ Ca2+ Mg2+

3.1655 2.584 2.361 1.398

0.1554 0.100 0.45 0.875

a Units: σ in Å and  in kcal/mol. Parameters of Na+ are from ref 35; parameters of Ca2+ and Mg2+ are from ref 34.

the ion was placed in the opposite direction at radial distance from 8 to 15 Å, providing ion-macroion separation from 11 to 18 Å, respectively. The dm shift of the macroion provides the minimal distance between the ion and the RRF boundary of the spherical cavity. This distance should exceed 12 Å; otherwise the interaction between the ion and the boundary can affect results of simulations.25 The water molecules were initially placed within a spherical volume with a radius of R ) 25 Å according to experimental water density. The radius of the initial exclusion sphere around the macroion was 10.8 Å. The mean force acting on the ion was evaluated for fixed positions of the ion with a step size of h ) 0.125 Å, implying 57 simulations for each mean force curve. The equations of motion of the rigid molecules were integrated by using quaternions, and the velocity version of the Verlet algorithm with a time step of 1 fs was used. The simulation time was tMD ) 0.5 ns, in addition to preceding equilibration. The velocities were scaled by using the procedure by Berendsen et al.39 with a time constant of 0.1 ps. The desired temperature was 298 K, and the averaged temperatures obtained were between 299 and 300 K. All MD simulations were carried out by employing the MOLSIM package.40 The potential of mean force at separation r between centers of the ion and macroion was integrated according to

u(r) ) u(rmax ) -

∫rr

F(ion)(r′) dr′

(1)

δi(ai ) ) δi - 1000NA

qi2 4π02waterai

γ

(3)

where δi is measured in M-1, qi is the ion charge, and γ is a coefficient from the decrease of permittivity of pure water in a weak external field water(E) ) water(0)(1 - γE2). Then the reduced values δi(ai ) are connected to the parameters ai and i via the equation22

δi(ai ,i) ) 1000NA

4π  33water(water - i) (a ) 3 i 2water + i

(4)

Parameters i and ai cannot be determined independently from eqs 3 and 4. A comparison of a PMF from MD simulations with finite-difference calculations of the electrostatic energy (Figure 2) makes it possible to adjust i, as the short-range repulsion of the ion becomes sensitive to this parameter when the macroion/water dielectric boundary intercepts the ion dielectric sphere. The reduced value δi(ai ) eliminates the solvent dielectric saturation effects originating from the field of the ion outside the ion dielectric sphere ai (Appendix 1).22 Therefore, the successful fitting of the parameters i and ai would justify neglecting this effect for the model shown in Figure 2. The parameters i, ai , and am were adjusted in such a way that (i) experimental δi values were satisfied via eqs 3 and 4 and (ii) u(r) from eq 1 was best fitted by the energy calculations for the continuum electrostatic model. The value u(rmax) in eq 1 was taken from the continuum electrostatic calculations. In those calculations, the dependence of SPC/E(E) on the electric field E of the macroion was calculated according to the Booth theory41 with the value SPC/E(0) ) 71 as described in ref 28. The function Booth(E(r)) of the distance r from the center of the macroion was obtained via the equation

max

by using the trapezoidal rule. Here F(ion)(r) is the mean force acting on the ion, and u(rmax) is the interacting potential at the maximal separation between centers of the ion and macroion (18 Å). The averaged force acting on the macroion was not counted due to its higher statistical uncertainty in comparison with F(ion). The continuum electrostatic model shown in Figure 2 was applied to adjust parameters am, ai , and i by means of a comparison of calculated electrostatic energies of this system with the MD-simulated PMF. The PPM used in these electrostatic calculations describes an ion as a sphere of low dielectric permittivity of ∼4 Å in radius with a point charge in the center. The sphere fills the volume inside the second coordination shell around the ion.22 As discussed in ref 22, parameters ai and i of the PPM are connected to experimental values δi for a decrease of the electrolyte solution permittivity with growth of the ionic concentrations ci19

solution(c) ) water -

∑i δici

(2)

First, these positive quantities δi should be “reduced” to the radius ai of the dielectric sphere approximating the polarization deficiency of the ion hydration shell (SI units)22

E(r) )

Qm 4π0Booth(E(r))r2

(5)

after 10 iterations starting from the dielectric constant of bulk SPC/E water (Booth(0) ) SPC/E ) 71); Qm is the macroion’s charge. Results. The all-atom (Figure 1) and dielectric (Figure 2) models were compared for three different values of the macroion charge Qm ) eZm. Values Zm ) 0, -8, and -16 were applied. In the case of Zm ) -16, the 8 Å radius macroion appears to be surrounded by a vast zone of saturated water, which cannot be neglected in the continuum dielectric calculations. The (r) dependencies according to eq 5 are shown in Figure 3. As follows from these results, the dielectric saturation effects can be disregarded at |Zm| < 8. Results of adjustment of the am, ai , and i parameters and values δi(ai ) are given in Table 2 for Ca2+, Mg2+, and Na+. In Figures 4-7, calculated electrostatic energies for the dielectric model are compared with PMFs from MD simulations. In all curves the value Zieψ(0) m (r) is subtracted from the potential. Here ψ(0) m (r) is the electrostatic potential of the isolated macroion for the dielectric model, i.e.,

ψ(0) m (r) )

Zme 4π0SPC/Er

(6)

Interaction between Two Double Helices of DNA

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5267

TABLE 2: Fitted Parameters of Dielectric Model and Derived δ Parameters ion

am/Å

ai/Å

i

δi(ai)a/M-1

b/M-1 δexp i

 b -1 δexp i (ai ) /M

Na+ Ca2+ Mg2+

10.3 10.3 10.0

3.8 4.0 4.1

25 7 7

9.4 16.5 17.8

11 (10)c 23 (21) 24 (22)

9.5 (8.5) 17.4 (15.4) 18.5 (16.5)

a Calculated for w ) 78.5 b Experimental data from ref 20. The first value for estimation δCl ) 2 M-1, and the value in brackets for δCl ) 3 M-1 (refs 19 and 22). c The experimental value δNaCl ) 12.8(4) from ref 21 is used.

Figure 5. Same comparison as in Figure 4, but for Qm ) 0. The large filled circles are MD simulation results. The solid and dashed curves are the ∆U for the electrostatic model with am ) 10.0 Å and am ) 9.7 Å, respectively. The small filled and open circles are MD results for Qm ) -16e and Qm ) -8e. Other parameters are as in Figure 4.

Figure 3. Dependence of the water permittivity Booth(r) around the spherical macroion in the absence of ions. The solid curve corresponds to Qm ) -16e, the dashed line is a result for Qm ) -8e.

Figure 6. Same comparison as in Figure 4, but for the Mg2+ ion. Qm ) -16e. The large circles are MD simulation results. The curve is ∆U for the electrostatic model with am ) 10.0 Å, ai ) 4.1 Å, i ) 7. The dotted curve is MD results for Ca2+ and Qm ) -16e.

Figure 4. Comparison of the PMF from MD simulations and the electrostatic energy of the dielectric model for Ca2+ and a charged macroion. ∆U(r) ) U(r) - qiΨm(r), where Ψm(r) is the potential of the macroion for water permittivity shown in Figure 3. Filled circles are MD results for Qm ) -16e, open circles are MD results for Qm ) -8e. The solid and dashed curves are the electrostatic energies of the dielectric model shown in Figure 2 for Qm ) -16e and Qm ) -8e, respectively. Parameters of the model are am ) 10.3 Å, ai ) 4 Å, i ) 7.

when the dielectric saturation of water is neglected and

ψ(0) m (r) )

Zme 4π0

∫r∞ 

dr 2 Booth(r)r

(7)

if the picture shown in Figure 3 is taken into account. In other words, the electrostatic energy of the dielectric model in excess to the electrostatic potential of the isolated macroion ∆uel(r) ) e/2 (Ziψm(r) + Zmψi(r)) - Zieψ(0) m (r) is compared to the similar quantity from MD simulations: ∆uMD (r) ) u(r) -Zieψ(0) m (r). Here u(r) is obtained from eq 1, and u(rmax) in eq 1 is determined as ∆uMD(rmax) ) ∆uel(rmax). As follows from adjustment of parameters am, ai , and i to fit the ∆uel(r) and ∆uMD(r) curves (Figures 4-7), values

Figure 7. Same comparison as in Figure 4, but for the Na+ ion. The filled and open circles are MD simulation results for Qm ) -16e and Qm ) -8e, respectively. The solid and dashed curves are ∆U for the electrostatic model with Qm ) -16e and Qm ) -8e, respectively. The model parameters are am ) 10.3 Å, ai ) 3.8 Å, i ) 25. The dotted curve is the same result at Qm ) -8e and am ) 10.8 Å.

δi(ai ,i) calculated via eq 4 agree with δi(ai ) values derived from experimental δi values through eq 3 (Table 2). It should be noted that the short-range repulsion of a cation from the macroion/water interface is quite sensitive to all three adjustable parameters. In other words, the PPM parameters ai and i not only fit the experimental δi for bulk electrolyte22 but also approximate the penetration of the ion hydration shell into the first layer of water molecules around the macroion. The

5268 J. Phys. Chem. B, Vol. 111, No. 19, 2007 repulsion of “structure-making” cations Na+ from the first layer of water agrees with early results of MD simulations of the air/ solution interface.42 Fitted radius am of the dielectric boundary satisfies the following dependence: am (Ca2+, Zm ) 0) < am(Mg2+, Zm ) -16) < am(Ca2+, Zm < -8) < am(Na+, Zm < -8). The value of am is between 10.0 and 10.8 Å, mainly ∼10.3 Å. This is in a good agreement with any expectation, as the LJ potential of the macroion is centered at r ) 8 Å and the first peak of hydrogen distribution around the isolated charged macroion is at r ) 10.0 Å and the first oxygen peak is at r ) 10.9 Å. The relationship between values of am under different conditions can be explained by the different destruction of the second coordination shell of the ion near the macroion/water interface. Indeed, the distribution function of solvent exceeds 2 in the second coordination shell of Mg2+ and Ca2+.43 In the case of Zm ) 0, water molecules are not so strongly polarized in the inner shell around the macroion as in the case of the central electric field, and the ion might retain some part of its second shell at short ion-macroion distances. It shifts the approximating dielectric boundary am toward the macroion. The Mg2+ ion possesses the stiffer second coordination shell, but the macroion’s field and high packing of water in the inner layer destroy it. It is more pronounced for Ca2+ and even more so for Na+. Moreover, the first coordination shell of Na+ starts to be destroyed in the high field of the macroion at low r (Figure 7, Zm ) -16). The main conclusion from this comparative analysis of ∆uel(r) and ∆uMD(r) is that the PPM can approximate not only the polarization deficiency of the hydrated ion in an external electric field but also the short-range repulsion of structuremaking ions from the water/nonpolarizable medium interface. Of course, this approximate model only imitates the true solventdriven force acting on an ion close to the interface. But the realistic position of the hard repulsion and relatively low magnitude of the solvent-driven oscillations in the simulated PMFs allow us to apply the model of a hydrated cation as a sphere of low permittivity to the development of the doublelayer theory based on the PPM. This will be done in the next section. Verification of the Continuum Dielectric Approximation of Solvent. In addition to the verification of the dielectric model of ion-macroion interactions, the continuum dielectric description of pure solvent around the charged spherical macroion should also be verified. To this aim, the electrostatic potential from eqs 6 and 7 was compared with the MD-simulated mean electrostatic potential at the oxygen atom of a SPC/E water molecule. The MD simulation model was as shown in Figure 1, but the macroion was fixed in the center of the spherical volume, and there was no ion in the system. Results for Zm ) -16 are shown in Figure 8 where the electrostatic potential ψ(0) m (r) for the dielectric model is compared to the mean potential from MD simulations: ψMD m (r) ) 〈ψw〉(r) + Zme/4π0r + (1/SPC/E -1)Zme/4π0RRF. Here 〈ψw〉(r) is the averaged electrostatic potential at the oxygen atom of a water molecule originating from the other molecules, and the last two terms are the electrostatic potentials of the macroion and its image charge, respectively. The thickness 2∆r of spherical layers to average 〈ψw〉(r) within r ( ∆r was 0.5 Å. After reaching the equilibrium, the time of simulation was 0.04 ns for each of 49 independent runs with different initial configurations. The final MD results appear as averaging over these 49 separate MD simulations.

Gavryushov

Figure 8. Comparison of the MD-simulated mean electrostatic potential acting on a water molecule with the continuum potential around the spherical macroion. Qm ) -16e. Filled circles are results of MD simulations. The solid curve is the potential for Booth(r) shown in Figure 3, the dashed curve is the potential of the macroion at (r) ) SPC/E. The thick dotted curve is the O-atom distribution and the thin dotted curve is the H-atom distribution around the charged macroion (the right scale).

As seen from Figure 8, the mean electrostatic force acting on a water molecule becomes oscillating mainly within the first coordination shell of the macroion. This agrees well with the reported results for the plain water surface.44,45 The disturbance of the “continuum” potential (eqs 6 and 7) is notable only within the inner solvent shell around the highly charged macromolecule. For the force acting on a charged ion, this oscillation effect will be diminished because the “frozen” hydration shell of the ion is retained even in high external electric fields and the diameter of the shell is ∼8 Å (Table 2). As a result, the mean electrical force acting on the ion from free water molecules is some average of the MD-simulated water potential (Figure 8) over an interval of a few a˚ngstroms. It is clearly seen from Figure 8 that such an average should be close to the “continuum” approximation (solid and dashed curves in Figure 8). This conclusion can explain the success of the continuum dielectric approximation (Figure 2) in application to the ion-macroion interaction. Despite the pronounced drop of the solvent permittivity around the highly charged macroion (Figure 3), the analysis of results shown in Figure 8 does not allow evaluation of this effect close to the macroion/water interface. And yet, this effect, confirmed by simulations of bulk water,45,22 cannot be eliminated for the spherical macroion as follows from fitting parameters ai and i in the case of Na+ and Zm ) -16 (Figure 7). The value of ∆uel(r) is always positive and asymptotically decreasing, as it is in the nature the charge-image repulsive interaction. After subtraction of ψ(0) m (r| ) SPC/E) (eq 6) from the MDsimulated PMF, ∆uMD(r) becomes asymptotically increasing, which exceeds the statistical errors (data not shown). It is eliminated by the use of eq 7 for ψ(0) m (r) (Figure 7). This effect is clearly observed for Na+ and Zm ) -16 but cannot be seen for divalent cations because the magnitude of ∆uel(r) is higher in this case. Epsilon-Modified Poisson-Boltzmann Equations In this section the theoretical foundations of a double-layer theory based on the PPM will be considered. The model of an ion of the double layer is shown in Figure 2. All assumptions, closures, etc. will be discussed in detail to derive the final set of equations. Model. As was shown in the previous section, a polarization model of the ion hydration shell can approximate the short-

Interaction between Two Double Helices of DNA

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5269

TABLE 3: E-MPB Parameters of Counterions aHC/Å

ion +

Na Ca2+ Mg2+

4.5 6.0 6.0

a

a/Å

i

3.8 4.0 4.1

25 7 7

a Estimated as a “mean” aHC in the approximate polarization model for NaCl (ref 26).

range repulsion of an ion from the interface between water and nonpolarizable medium. Since the model also leads to a good approximation of the ion energy in the external electric field,22 this description of ions can be used in the statistical mechanical equations to explore polarization effects in a dense double layer. An attempt to incorporate such effects in the MPB equations was done in ref 28, where it was demonstrated that the effects could be significant in the case of DNA electrostatics and 2:1 electrolyte. However, the equations applied in that work were only approximate ones, and it will be shown below that they are correct only for weak fields and dilute electrolyte solutions. Recent improvements of the adjustable parameters of ion hydration shells and ion-ion interactions (refs 22, 25, and 26; Table 2) allow us to apply a more sophisticated theory to systems of interest. The PPM parameters of hydrated Na+, Ca2+, and Mg2+ are collected in Table 3 from results of the previous section and refs 25 and 26. We simplify the theory by the case when only cations are counterions in the double layer. Then properties of the structure-breaking42 co-ions Cl- can be disregarded, and low polarization deficiency of the Cl- hydration shell19,22 can be neglected. The ion-ion hard-core interaction distance aHC ij is set to a common value of aHC that is correct for divalent cations or chosen as a mean for both Na+ and Cl-.26 It will be important for the following that aHC > ai for all alkali and alkaline earth metal cations. There is another important inequality related to distribution function gi(r) of ions of species i. Since only one ion center can be found within the sphere |r′ - r| < aHC/2, the upper limit of gi(r) is determined by the expression

∑i ni0 ∫|r′-r| ai of solvent around the ion (Appendix 1). Therefore, any difference between double-layer properties predicted for nonpolarizable ions and for ions shown in Figure 2 could be more pronounced in real systems. “Mean” Electrostatic Potential. To derive a closed set of equations for a cloud of such ions as shown in Figure 2, we will have to define some quantities. Let the thermodynamic values be averaged when an ion of species i appears at point r1. Then the mean electrostatic potential at point r will be ψi(r1,r). This potential can be divided into several components F

ψ(r1,r) ) φ (r1,r) + φ (r1,r) + φc (r1,r) + ψ*(r1,r) δ



(9)

(In the following the ion subscript will be omitted in some expressions.) Here φδ(r1,r) ) ψ(r1,r|qi) - ψ(r1,r|0) is the potential originating from the presence of the central charge of the ion

(parameters ai and i are treated as disengaged), φ(r1,r) ) ψ(r1,r|0,i) - ψ(r1,r|0,(r)) is the disturbance of the potential caused by polarization of the medium due to the ion at point r1 ((r) denotes the mean permittivity of the solution as it will be defined below), φcF(r1,r) ) ψ(r1,r|0,(r)) - ψ*(r1,r) is a potential change caused by the disturbance of the mean ionic distributions at |r - r1| > aHC due to the hard-core sphere at point r1, and ψ*(r1,r) is the rest of the potential. According to these definitions, ψ*(r1,r1) is a potential in the center of sphere |r - r1| < aHC, and it averages potentials originating from all free charges outside this sphere at undisturbed mean ionic number densities nk(r) ) nk0gk(r). Let φk(r,r1) denote the mean potential at point r1 created by an ion of species k fixed at point r. (It should be distinguished from ψk(r,r1) that is the mean potential at point r1 from all charges when ion k is found at point r.) Then

ψ*(r1,r1) )

∑k ∫|r-r |>a 1

HC

φk (r,r1)nk(r) dr

(10.1)

or, denoting φk(r,r1) ) qkGk(r,r1)

ψ*(r1,r1) )

∑k ∫|r-r |>a

HC

1

Gk(r,r1)nk(r)qk dr (10.2)

For function φk(rq,r) at |r - rq| > aHC, one can write the following relations for the mean electric field of ion k Ek(rq,r) ) -∇φk(rq,r) and electrical displacement Dk(rq,r)

Ek(rq,r) )

∑j ∫|r-r′| aHC 0, outside, |r - r1|> aHC, |r - r′|< aHC - qiδ3 (r - r1)/0, outside, |r - r1|< aHC - Fm (r)/0, inside (21)

For the dependence (r1,r) we assume that

(r1,r) ≈

{

i, |r - r1|< a (r), |r - r1|> a

(22)

As for eq 13, it is convenient to introduce potentials similar to ψ ˜ (r) and φF0(r′,r) from eq 14: ψ ˜ i (r1, r) ) ψ ˜ /i (r1,r′,r) F / ˜ i (r1, r, r) - ψ ˜ *(r,r) from eq φi0(r1,r′,r). Then the difference ψ ˜ *(r,r) ) ψ ˜ i(r1,r) - ψ ˜ (r) 20.2 will be written: ψ ˜ /i (r1,r,r) - ψ F (r1,r,r) - φF0(r,r)). Further we admit that at |r - r1| > + (φi0 aHC F ˜ (r) | . | φi0 (r1,r,r) - φF0(r,r)| |ψ ˜ i(r1, r) - ψ

(23)

Interaction between Two Double Helices of DNA

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5271

which yields the important expression

Wik (r1,r) - Wk(r) ≈ uik (|r - r1|) + ˜ i (r1,r) - ψ ˜ (r)) (24) qk (ψ The boundary value problem for ψ ˜ i(r1,r) is

˜ i(r1,r)) ) ∇((r1,r)∇ψ

{

∑

nik(r1,r)qk/0, outside, |r - r1|> a , -qiδ3(r - r1)/0, outside, |r - r1|< aHC, (25) -Fm (r)/0, inside HC

where (r1,r) satisfies eq 22. By subtracting eq 15 from eq 25 and using eqs 22 and 24, we obtain the boundary value problem for the fluctuation potential φi(r1,r) ) ψ ˜ i(r1,r) -ψ ˜ (r)

0∇((r)∇φi(r1,r)) ) -

∑k

(nik(r1,r) - nk (r))qk )

(exp(-β(Wik(r1,r) - Wk(r))) - 1) ≈

|r - r1|> aHC or inside ∇((r1,r)∇φδi (r1,r)) ) -qiδ3(r - r1)/0,

(29.1)

|r - r1|< aHC, outside

On introducing φi0(r1,r) ) φi(r1,r) + φiF(r1,r), we have the second problem from eqs 26 and 28

∇((r)∇φi0(r1,r)) - κ2(r)φi0(r1,r) ) 0, |r - r1|> aHC or inside (29.2) ˜ (r)), ∇((r1,r)∇φi0 (r1,r)) ) -∇((r1,r)∇ψ |r - r1|< aHC, outside

nk(r) ) nk0 exp(-βWk(r))

∑k nk(r)qk2 βφi(r1,r)

|r - r1| > aHC, outside the macromolecule; β ) 1/kT; ∇((r1,r)∇φi(r1,r)) ) -∇((r1,r)∇ψ ˜ (r)) - qiδ3(r - r1)/0, HC |r - r1| < a , outside the macromolecule; ∇(m∇φi(r1,r)) ) 0 inside the macromolecule. The expression for |r - r1| > aHC is the linearized Loeb’s closure as used in the MPB theory for nonpolarized ions,47 and the expression for |r - r1| < aHC is the definition of φi(r1,r) substituted in eq 25. Finally, the boundary value problem for the fluctuation potential can be written as

{

{ {

∇((r)∇φδi (r1,r)) - κ2(r)φδi (r1,r) ) 0,

where ψ ˜ (r) satisfies eq 15, (r1,r) is described by eq 22, (r) is connected to nk(r) via eq 19, and nk(r) is connected to the ion PMF as

∑k nk0(exp(-βWik (r1,r)) - exp(-βWk(r)))qk )

∑k nk0qk exp(-βWk(r)) ×

From the definition of φδi (r1,r) (eq 9) and eqs 26 and 28, it is evident that φδi (r1,r) is a solution of the boundary value problem

∇((r)∇φi(r1,r)) - κ2(r)φi(r1,r) ) 0, |r - r1|> aHC or inside ∇((r1,r)∇ φi(r1,r)) ) -∇((r1,r)∇ ψ ˜ (r)) - qiδ3 (r -r1)/0, (26) |r - r |< aHC, outside 1

where κ2(r) ) β ∑nk(r)qk2/0. According to definitions, ψi(r1,r1) ) ψ/i (r1,r1,r1) ) / ˜ i(r1,r1), therefore ψ ˜ i (r1,r1,r1) ) ψ

ψi(r1,r1) ) ψ ˜ i(r1,r1) ) φi(r1,r1) + ψ ˜ (r1)

(27)

Omitting index i and using the definition of ψ ˜ (r), we obtain ψ(r1,r1) ) φ(r1,r1) + ψ*(r1, r1) - φF0(r1,r1). Equation 9 written for r ) r1 is ψ(r1, r1) ) φδ (r1,r1) + φ(r1,r1) + φFc (r1,r1) + ψ*(r1,r1). After the substitution of ψ(r1,r1) from the previous expression, it leads to the equation φ(r1,r1) ) φδ(r1,r1) + φ(r1,r1) + φFc (r1,r1) + φF0(r1,r1). Denoting φF(r1,r1) ) φFc (r1,r1) + φF0(r1, r1) and writing ion index i again, we have the final expression

φi(r1,r1) ) φδi (r1,r1) + φi (r1,r1) + φFi (r1,r1)

(28)

(30)

After establishing the final link between Wi(r) and the potential ψi(r,r), the set of eqs 15, 19, 22, 27-29, and 30 will be closed, and it can be solved by an iterative process as described in refs 11 and 28. This link follows from the Kirkwood charging process.47 The ion PMF can be divided into the two parts Wi(r) ) W0i (r) + Wqi (r). The first term is a work of transfer of an uncharged ion (with all its dielectric properties) from infinity to the point r, and the second term is a work of charging the ion at the point r. The term W0i (r) includes the work against the osmotic pressure of ionic gas and a change of the electrostatic potential of the macromolecule’s charges caused by the spherical cavity in the double layer and by the ion dielectric sphere. From eqs 27 and 28 and the definition of φ0i (r1,r), the mean electrostatic potential in the center of the ion at point r1 is written as ψi(r1,r1) ) φδi (r1,r1) + φ0i (r1,r1) + ψ ˜ (r1), so the PMF Wi(r1) is

Wi(r1) ≈ β-1

(

∑k ∫|r-r | 8.8 Å) was set to the bulk value depending on the bulk ionic concentrations (eq 2) and the ion exclusion zone started at r < 10.5 Å.48 Comparisons of results of -MPB and ordinary MPB calculations are shown in Figures 9-11. For 2:1 electrolyte, both MPB and -MPB results are notably different from those of the PB calculations (Figure 10) due to the finite size of ions and interionic correlation effects.47,11 In the case of 1:1 electrolyte (NaCl), the PB equation is a good approximation (Figure 11). In this case magnitudes of both W0i (r) and WiCSI(r) ) 0.5 q2i (φδi (r) - φδi (∞)) + qiφ0i (r) do not includes exceed 1 kT (data not shown). (The PMF term WCSI i

Interaction between Two Double Helices of DNA

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5273 Ion Polarization Effects on Energy of DNA-DNA interactions

Figure 10. Mean electrostatic potential around the cylindrical model of DNA in CaCl2 solutions. Thick curves are obtained for high concentration of electrolyte (c2 ) 0.167 M) and thin curves for c1 ) 0.0167 M. “1D” solutions. Solid curves the -MPB equations, dashed curves the MPB equations (nonpolarizable ions), dotted curves the PB equation.

To explore the possible effects of the ion hydration shell polarization deficiency on the interactions between two macromolecules of DNA, the 3D -MPB equations were solved for a system consisting of two parallel cylinders at the axis separation d. The parameters of the models of DNA and ions are described in the previous section. The free energy of the system of two cylinders was calculated for different d varying from 19 to 40 Å and for the bulk concentrations c1 and c2 of NaCl and CaCl2. Results were compared with those of the MPB calculations for nonpolarizable ions and of the PB equation. The free energy was obtained via two processes: a transfer of the uncharged macromolecules into electrolyte solution and charging them up28

W ) W0 + Wch

(32)

where

Wch )

Figure 11. Same as in Figure 10, but for NaCl. The high concentration of electrolyte is c2 ) 0.5 M and the low c1 ) 0.05 M. Notation as in Figure 10.

Figure 12. Drop of the local permittivity (r) around the cylindrical model of DNA. “1D” solutions of the -MPB equations. Thick curves are results for c2 concentrations of CaCl2 and NaCl, and thin curves for c1. Solid curves are for CaCl2, and dashed curves for NaCl.

effects of interionic correlations, ion finite size, and ion-image interaction.)28 The drop of permittivity close to the cylinder is shown in Figure 12. Despite the improved parameters of the dielectric model of ions, this picture is similar to the previously published -MPB results.28,22 Results of the “1D” and full “3D” numerical solutions of the -MPB equations are compared in Figure 9 for c1. They are almost identical, which is also observed for c2 results (data not shown). It allows us to apply the 3D -MPB code to a study of DNA-DNA interactions in the solutions of CaCl2 and NaCl.

∫01 dλ ∫ dr Fm(r)ψ(r,λ)

(33)

and W0 ) β-1 ∑ ∫ (n0k - nk (r)) dr. Here ψ(r,λ) is the mean potential from a solution of the MPB equations at charge density of the macromolecule equal to λFm(r). The energy calculation was intensively tested for the one cylinder model described in the previous section. By a comparison of “3D” and “1D” solutions obtained with different number of charging states, it was shown that only 3 to 4 charged states provided a good accuracy of eq 33. The model of the two cylinders requires some more comments. The solution permittivity (r) obeying eq 19 starts to separate the two cylinder interiors at d > 2 × am ) 2 × 8.8 Å ) 17.6 Å. It looks reasonable because the minimal distance between the outer atom centers of the cylinders is 17.6 - 2 × 6.5 ) 4.6 Å that is a bit less than the doubled diameter of a water molecule (2 × dw ) 2 × 2.8 Å ) 5.6 Å). The ion center can pass between the cylinders at d > 2 × 10.5 Å ) 21 Å, i.e., the “free space” 21 - 2 × (6.5 + dw/2) ) 5.2 Å that is a bit less than about 6.5 Å diameter of the ion with water molecules of its inner shell of hydration.43 So, this model is expected to describe ions and water around the two cylinders as well as the dielectric approximation fits results of MD simulations in Figures 4-7. The DNA model is apparently a very rough approximation, but in this study we estimate only the ion polarization effects, i.e., the difference between -MPB and MPB results. Results of -MPB, MPB and PB calculations of the free energy G12 for the model of the two cylinders are shown in Figures 13-14 as ∆G(d) ) G12(d) - G12(∞). Results for the bulk concentration c1 (0.05 M of NaCl and 0.0167 M of CaCl2) are collected in Figure 13, whereas results for c2 (0.5 and 0.167 M, respectively) are shown in Figure 14. For MPB calculations with nonpolarizable ions, an attraction of the cylinders is observed at d > 27 Å for both concentrations of CaCl2. This result is in a good agreement with reported MC simulations for the PM ions where electrolyte-induced attraction of polyanions of DNA was observed.49-51 This attraction can be explained by interionic correlations within the double layer. The attraction of the cylinders disappears in the case of the polarizable hydration shells of ions and the -MPB calculations lead to results that are closer to ones from the PB equation. In this case another feature is observed: the solution of the -MPB

5274 J. Phys. Chem. B, Vol. 111, No. 19, 2007

Figure 13. Free energy of electrostatic interaction between two cylindrical models of DNA at low bulk concentrations c1 of NaCl and CaCl2. Circles are results for 0.0167 M CaCl2, and triangles for 0.05 M NaCl. The filled symbols are results of -MPB calculations; the open symbols are results of MPB calculations for nonpolarizable ions. The solid curve is results of PB calculations for CaCl2, and the dashed curve PB results for NaCl.

Figure 14. Same as in Figure 13 but for c2 (0.167 M CaCl2 and 0.5 M NaCl). Notations as in Figure 13.

equations does not exist at d < 24 Å for CaCl2 and at d < 21 Å for NaCl. This requires some comments. The less d the more iterations of the main loop of eqs 15-29-31-15 are needed to obtain a convergence of the -MPB solution for the fully charged state of the two cylinders. At d ) 25 Å, this number reaches 16 for the c2 CaCl2 calculations. The final drop of (r) is shown in Figure 15. At d < 24 Å, the function (r) between the two cylinders ceases to converge and decreases until reaching the lower value i. Even before that, the inequality in eq 8 becomes broken in this area, i.e., we cannot consider ions as freely moving in some area between the cylinders. In other words, at d < 24 Å ions together with their distorted hydration shells become bound in this central zone. This result does not contradict with the phenomenological theory of short-range hydration forces between macromolecules and bilayers.52-54 One of the important aspects of the hydration forces is the role of ions.55,56 In the phenomenological theory,54 the role of ions is reduced to their adsorption within the inner core around the macromolecule, whereas the mean electrostatic potential obeys the linear PB equation outside. The nonexistence of -MPB solutions at d < 24 Å (Figures 13 and 14) predicts such adsorption of Ca2+ in the strong field of two DNA macromolecules. It should be noted that this result was obtained for a simplified cylindrical model, whereas for an all-atom geometry model of DNA, both MPB calculations11 and MC simulations for PM ions51,57 suggest that the ionic distributions are sensitive to real geometry of the macromolecule.

Gavryushov

Figure 15. Contour plot of (r) for the model of two cylinders of DNA in 0.167 M CaCl2 solution. The interaxial distance d ) 25 Å. The ticks of the plot boundary are separated by 5 Å. The contours of  are 70, 50, 30, and 10 from the outer zone to the interior of the two cylinders. The inner dielectric constant m of the 8.8 Å radius core of the cylinders is 2. The bulk permittivity of the electrolyte solution (c2) ) 75.8.

On the other hand, both Na+ and Ca2+ do not prevent indefinite swelling of DNA,56 which is confirmed by the present -MPB calculations. Another experimental fact found for MgCl255 is also clearly observed in Figures 13 and 14 for CaCl2: the shorter decay length λ than λDH from the PB calculations and almost independence of the DNA-DNA force on the electrolyte concentration. Similar experimental phenomenon for NaCl55 are not so much pronounced from the -MPB results, but it should be mentioned that the final conclusion could only be drawn from -MPB calculations for a hexagonal packing of all-atom geometry models of DNA, as the cell model55 cannot be applied in MPB calculations at low d due to the fluctuation potential problem (eqs 29). Sufficient concentrations of some transition-metal ions and trivalent ions cause spontaneous DNA assembling into ordered arrays of helices with a surface separation of 7-12 Å depending on counterion species but not its concentration over the critical level.58,59 Such an attraction is observed in MPB results for nonpolarizable divalent counterions (Figures 13 and 14). However, it should be regarded rather as a coincidence because of the DNA model simplifications. For the same reason, the -MPB results for Ca2+ shown in Figures 13 and 14 should not be regarded as a “prediction” of DNA swelling in the presence of Ca2+ ions either. The fact of assembling DNA in the presence of Cd2+ or Mn2+ ions at higher temperatures56,59 might suggest either (i) notably different polarization properties of their hydration shells than those of Ca2+ (and easier destruction of the hydration shells of Cd2+ and Mn2+ at higher temperatures) or (ii) an importance of ionic adsorption in the grooves of real DNAsthe effect naturally expected to be ion-specific. The 3D -MPB calculations for all-atom geometry model of DNA could verify the latter assumption, although more reliable results should only be obtained by all-atom MD simulations of DNA with ions. The present comparisons of -MPB with MPB calculations just indicate that polarization of hydration shells of divalent cations can be a notable contributor in DNA-DNA interactions. Conclusions As follows from comparisons with all-atom MD simulations reported in the present work, the PPM of ions well approximates the short-range repulsion of cations Na+, Ca2+, and Mg2+ from the interface between water and a nonpolarizable medium at different intensities of the applied external field. Besides, the fitted PPM parameters of ions lead to predictions of the

Interaction between Two Double Helices of DNA experimental dielectric constants of bulk electrolyte solutions. This allows us to apply the PPM to a double-layer theory in order to explore the role of solvent dielectric saturation in electrostatics of macromolecules. The very idea to incorporate polarization deficiency of the ion hydration shell into the MPB theory of double layers was discussed nearly a half a century ago.29-32 An attempt to apply such -MPB equations to DNA electrostatics was done in ref 28. In the present work, the -MPB equations are derived on a deeper theoretical basis and all assumptions and closures are considered in details. In addition to the linearized Loeb’s closure and the Kirkwood hierarchy of distribution functions, the -modified equations are based on approximations expressed by eqs 18, 20.1, 22 and 23. As follows from eq 29.2, the field of the polarization charges of the ion hydration shell is screened like any electrostatic interaction in the cloud of counterions. On the contrary, the dependence of the local solution permittivity (r) on the ionic concentration remains to be unscreened (eq 19).60 In other words, these two effects are different at finite ionic concentrations, although they are physically undistinguishable at infinity dilution of electrolyte. Outside the macromolecule, the “mean” potential ψ(r) from the Poisson eq 15 is not the true mean potential but only its part coming from outside the ion hardcore repulsion sphere centered at r. Applications of the improved -MPB equations to simplified models of one and two DNA macromolecules have shown that (i) the ion hydration shell polarization effects cannot be neglected even for diluted solutions of 2:1 electrolyte and (ii) at high intensities of the electric field, the -MPB calculations predict the bound state of ions in the inner shell around the macromolecule. Some of -MPB results for the two simplified macromolecules of DNA in solution of CaCl2 bear close resemblance with experimental data of DNA-DNA force reported by Rau and co-workers.56 The further verification of the theory should be obtained from -MPB calculations for allatom geometry model of DNA and macromolecular configuration matching exactly conditions of the experiments. Appendix 1 Polarization Energy of an Ion Hydration Shell outside an Ion Dielectric Sphere. The outer part of the polarization energy RiE2 of the ion hydration shell in the external field E can be written as

∆uipol,outer(E) )

(

qi2

)

2

qi 1 1 ≈ γE2  (E)  (0) 8π0ai 8π0(0)ai (A.1)

The latter expression can be obtained not only by the ion charging process as written in eq A.1 but also by the increase E of the electric field E.22 Then ∆upol i (E) ) -∫0 ∆Pi(E′) 3(E′)/ (2(E′) + 1) dE′, where ∆Pi(E) is the difference between the dipole moment Pi(E) of the infinitely large volume of water with one ion “i” and the dipole moment P(E) of the same volume of pure water. At the outer volume integration, it should be taken into account that the projection of polarization pi(E,r, θ) onto direction of E will be pi(E,r,θ) ) 0((E,r,θ)-1)Ei(E,r, θ) and

(E) ) lim 2π ∆Pouter i Rf∞

∫aR dr r2 ∫0π dθ sin θ ×  i

(pi (E,r,θ) - 0 ((0)(1 - γE2) - 1)E)

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5275 where

(

(E,r,θ) ≈ (0) 1 - γE2 - 2Eγ

qi 4π0(0)(1 - γE2)r2

)

cos θ

and

Ei(E,r,θ) ≈ E +

qi 4π0(0)(1 - γE2)r2

cos θ

Here θ is the angle between vectors E and r. For the bulk solution symmetry, only the quadratic term proportional to q2i cos2 θ leads to a nonzero term in ∆Pi(E), which results in final (E) originates from expression (A.1). Thus, this part ∆upol,outer i nonlinear polarization of the solvent in the sum of electric fields E and that of the ion. The terms linear in respect of cos θ will not be cancelled if the ion is located close to the water/ macromolecule interface or within a non-homogeneous external field. They are, however, expected to be much lower than a value in eq A.1 because they are a change of the image interaction energy of the ion charging process at replacement of (0) by (E). Anyway, one can conclude from eq 3 (or, equally, from eq A.1) and from Table 2 that in the symmetrical case of bulk electrolyte solution, the outer part δi - δi(ai ) is much lower than δi(ai ) originating from the inner saturated zone of the ion hydration shell (eq 4). This “outer” part of the ion free energy will be neglected for the model shown in Figure 2, so the PPM actually describes less polarizable hydration shells of ions than real ones. Acknowledgment. The author thanks Per Linse for the MOLSIM program and for discussions. Computation resources of the Moscow Joint Supercomputer Center are gratefully acknowledged. References and Notes (1) Cheatham, T. E., III.; Kollman, P. A. J. Mol. Biol. 1996, 259, 434. (2) Young, M. A.; Jayaram, B.; Beveridge, D. L. J. Am. Chem. Soc. 1997, 119, 59. (3) Young, M. A.; Ravishanker, G.; Beveridge, D. L. Biophys. J. 1997, 73, 2313. (4) Sprous, D.; Young, M. A.; Beveridge, D. L. J. Phys. Chem. B 1998, 102, 4658. (5) Long, H.; Kudlay, A.; Schatz, G. C. J. Phys. Chem. B 2006, 110, 2918. (6) Gil Montoro, J. C.; Abascal, J. L. F. J. Chem. Phys. 1995, 103, 8273. (7) Lyubartsev, A. P.; Laaksonen, A. J. Chem. Phys. 1999, 111, 11207. (8) Fixman, M. J. Chem. Phys. 1979, 70, 4995. (9) Jayaram, B.; Sharp, K. A., Honig, B. Biopolymers 1989, 28, 975. (10) Das, T.; Bratko, D.; Bhuiyan, L. B.; Outhwaite, C. W. J. Chem. Phys. 1997, 107, 9197. (11) Gavryushov, S.; Zielenkiewicz, P. Biohys. J. 1998, 75, 2732. (12) Bacquet, R. J.; Rossky, P. J. J. Phys. Chem. 1984, 88, 2660. (13) Gonzales-Tovar, E.; Lozada-Cassou, M.; Henderson, D. J. Chem. Phys. 1985, 83, 361. (14) Boda, D.; Gillespie, D.; Nonner, W.; Henderson, D.; Eisenberg, B. Phys. ReV. E 2004, 69, 46702. (15) Skolnick, J.; Fixman, M. Macromolecules 1978, 11, 867. (16) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219. (17) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415. (18) Dong, F.; Vijayakumar, M.; Zhou, H.-X. Biophys. J. 2003, 85, 49. (19) Hasted, J. B.; Ritson, D. M.; Collie, C. H. J. Chem. Phys. 1948, 16, 1. (20) Harris, F. E.; O’Konski, C. T. J. Phys. Chem. 1957, 61, 310. (21) Buchner, R.; Hefter, G. T.; May, P. M. J. Phys. Chem. A 1999, 103, 1. (22) Gavryushov, S.; Linse, P. J. Phys. Chem. B 2003, 107, 7135. (23) Levine, S.; and Rosenthal, D. K. In Chemical Physics of Ionic Solutions; Conway, B. E., Barradas, R. G., Eds.; Wiley: New York, 1966.

5276 J. Phys. Chem. B, Vol. 111, No. 19, 2007 (24) Levine, S.; Bell, G. M. In Electrolyte Solutions; Pesce, M., Ed.; Pergamon: New York, 1962. (25) Gavryushov, S.; Linse, P. J. Phys. Chem. B 2006, 110, 10878. (26) (a) Gavryushov, S. J. Phys. Chem. B 2006, 110, 10888. (b) Gavryushov, S. J. Phys. Chem. B 2006, 110, 13678. (27) Lamm, G.; Pack, G. R. J. Phys. Chem. B 1997, 101, 959. (28) Gavryushov, S.; Zielenkiewicz, P. J. Phys. Chem. B 1999, 103, 5860. (29) Prigogine, I.; Mazur, P.; Defay, R. J. Chim. Phys. 1953, 50, 146. (30) Bolt, G. H. J. Colloid Sci. 1955, 10, 206. (31) Sparnaay, M. J. Recl. TraV. Chim. Pays-Bas. 1958, 77, 872. (32) Levine, S.; Wrigley, H. E. Discuss. Faraday Soc. 1957, 24, 43. (33) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (34) Åqvist, J. J. Phys. Chem. 1990, 94, 8021. (35) Dang, L. X. J. Am. Chem. Soc. 1995, 117, 6954. (36) Linse, P. J. Phys. Chem. 1986, 90, 6821. (37) Rami Reddy, M.; Berkowitz, M. Chem. Phys. Lett. 1989, 155, 173. (38) Svishchev, I. M.; Kusalik, P. G. J. Phys. Chem. 1994, 98, 728. (39) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (40) Linse, P. MOLSIM, version 3.0; Lund University: Lund, Sweden 2000. (41) (a) Booth, F. J. Chem. Phys. 1951, 19, 391. (b) Booth, F. J. Chem. Phys. 1951, 19, 1327. (c) Booth, F. J. Chem. Phys. 1951, 19, 1615. (42) (a) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2001, 105, 10468. (b) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361. (c) Mucha, M.; Frigato, T.; Levering, L. M.; Allen, H. C.; Tobias, D. J.; Dang, L. X.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 7617. (d) Manciu, M.; Ruckenstein, E. Langmuir 2005, 21, 11312. (43) Babu, C. S.; Lim, C. J. Phys. Chem. B 1999, 103, 7958.

Gavryushov (44) Manciu, M.; Ruckenstein, E. Langmuir 2005, 21, 11749. (45) Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 110, 7935. (46) Haggis, G. H.; Hasted, J. B.; Buchanan, T. J. J. Chem. Phys. 1952, 20, 1452. (47) Carnie, S. L.; Torrie, G. M. AdV. Chem. Phys. 1984, 56, 141. (48) The value of 10.5 Å follows from MD results in Figures 4-7. In that case the hard repulsion of ions starts at r ≈ 12 Å., but the outer atom centers of the macroion are at am ) 8 Å. In the case of the DNA cylinder, am ) 6.5 Å, leading to estimation of 10.5 Å. (49) Guldbrand, L.; Nilsson, L. G.; Nordenskio¨ld, L. J. Chem. Phys. 1986, 85, 6686. (50) Allahyarov, E.; Lo¨wen, H. Phys. ReV. E 2000, 62, 5542. (51) Allahyarov, E.; Gompper, G.; Lo¨wen, H. Phys. ReV. E 2004, 69, 41904. (52) Marcˇelja, S.; Radic´, N. Chem. Phys. Lett. 1976, 42, 129. (53) Kornyshev, A. A.; Leikin, S. Phys. ReV. A 1989, 40, 6431. (54) Kornyshev, A. A.; Leikin, S. J. Chem. Phys. 1997, 107, 3656. (55) Rau, D. C.; Lee, B.; Parsegian, V. A. Proc. Natl. Acad. Sci. U.S.A. 1984, 81, 2621. (56) Leikin, S.; Parsegian, V. A.: Rau, D. C.; Rand, R. P. Annu. ReV. Phys. Chem. 1993, 44, 369. (57) Abascal, J. L. F.; Gil Montoro, J. C. J. Chem. Phys. 2001, 114, 4277. (58) Rau, D. C.; Parsegian, V. A. Biophys. J. 1992, 61, 246. (59) Rau, D. C.; Parsegian, V. A. Biophys. J. 1992, 61, 260. (60) This phenomenon was observed in MC simulations of CaCl2 and MgCl2 activity coefficients with ion-ion potentials obtained from all-atom MD simulations. The best adjustment of the experimental data required screening of the long-range polarization repulsion of ions together with experimental dependence (c) from eq 2 (ref 25).