Monte carlo simulations of counterion accumulation near helical DNA

Monte carlo simulations of counterion accumulation near helical DNA

Volume 129, number 2 CHEMICAL PHYSICS LETTERS MONTE CARLO SIMULATIONS OF COUNTERION ACCUMULATION NEAR HELICAL DNA Pamela MILLS, Mark D. PAULSEN, ...

402KB Sizes 2 Downloads 23 Views

Volume 129, number 2

CHEMICAL PHYSICS LETTERS

MONTE CARLO SIMULATIONS OF COUNTERION ACCUMULATION

NEAR HELICAL

DNA

Pamela MILLS, Mark D. PAULSEN, Charles F. ANDERSON Departments

22 August 1986

and M. Thomas RECORD Jr.

of Chemistry and Biochemrstty, Unrversrty of Wisconsin at Madison, Maduon, WI 53706, USA

Received 1 April 1986; in final form 10 June 1986

For a model of DNA having a helical charge distribution, the number of adjacent counterions is calculated from the Monte Carlo (MC) simulations at various salt concentrations. These results are compared with corresponding MC simulations for an axial polyion charge distribution and with the predictions of other polyelectrolyte theories.

1. Introduction The Monte Carlo (MC) method has been widely used to compute ion distributions surrounding model polyions of various geometries [l-9]. These studies have been directed primarily at testing the molecular and thermodynamic predictions of approximate analytical polyelectrolyte theories, in particular the Poisson-Boltzmann (PB) theory. MC results for univalent ion distributions around cylindrical model poly ions having some of the average structural characteristics of native BDNA are generally close to the corresponding distributions calculated with the PB equation [2-41. Moreover, these MC simulations confirm the prediction based on PB calculations that the number of counterions near DNA (those, for example, found within 0.4 nm of the polyion surface) varies significantly with the total salt concentration over the typical experimental range. In contrast, for the same model of the polyelectrolyte solution the counterion condensation (CC) hypothesis predicts that the extent of counterion association with DNA is salt-invariant [lo,1 11. This prediction is in accord with a two-state interpretation of data from various NMR experiments [ 12-171. The simple axial model of the DNA charge distribution, which we investigated in paper I [4], is replaced in the work reported here by a more realistic helical charge distribution. Our objective is to determine whether this refinement of the polyion model can resolve the apparent conflict between experimen0 009.2614/86/$ 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

tal observations and previously reported MC predictions for the salt dependence of the number of counterions in the close vicinity of the DNA.

2. Model and method

The Metropolis Monte Carlo method, which we applied in paper I to a model polyion with an axial charge distribution, has been used here to simulate ion distributions around a model polyion with a helical charge array. Specifically, two full turns of the DNA helix were modeled as two dyadic helical strands (180’ out of phase), each containing 20 univalent charges, on the surface of an impenetrable cylinder of radius a0 = 1.Orun. The projected axial charge spacing is b = 0.17 run. For all the simulations reported here, the outer boundary of the cylindrical MC cell is R = 10 MI, which corresponds (in the usual cell model approximation) to a DNA phosphate concentration of 0.031 mol/dm3. The univalent counterions and coions were modeled as spheres of equal diameter 6, and the aqueous solvent was modeled as a uniform dielectric medium. Two values of 6 (0.1 and 0.6 run) were investigated to test for any effect of ion size on the distribution of counterions near the surface of the helical DNA. All simulations were performed for a temperature of 298 K, and a dielectric constant E of 78.5. The energy of a particular configuration of ions in 155

.the MC cell, calculated from Coulombic interparticle potentials, can be expressed as the sum of three types of terms: N

40

NN

where N is the total number of counterions and coions in the MC cell, Uzr is the pairwise interaction energy for the mobile ions (eq. (2) in paper I), @’ is defined analogously for the interactions between F” rxed polyion char es and mobile ions; and the external potential @Ts9 , (unlike ‘P,t in paper I) also has contributions from the portions of the polyion external to the MC cell. In calculating @ST, the polyion charges outside the MC cell are approximated as finite lines of uniform, continuous charge, and the contribution of the mobile ions is estimated during the simulation (with eqs. (5) and (6) of paper I) by using periodically updated cumulative histograms for the ion distributions inside the MC cell. In reality the small-ion distributions near DNA must exhibit some dependence on angular and axial coordinates, because the phosphates on DNA are arranged in a discrete helical array. However, previous MC [2] and PB [18] calculations of ion distributions surrounding a model polyion having a helical charge distribution indicate that they in effect are exclusively radial functions at distances greater than aO.5 mn from the polyion surface. Thus, in our simulations the average charge distribution in the central cell, and hence @zT, are approximated as purely radial functions. The construction of transition probabilities using AE, calculated with eq. (l), and all other aspects of our computations follow the standard Metropolis procedures used previously [3,4]. Long-range Coulombic interactions produce a high density of counterions in the volume of solution near a highly charged polyion, such as DNA. This purely electrostatic effect can be expressed in terms of p(A), the number of counterions per polyion charge within an annular volume of thickness A surrounding the polyion: a+A P(A) = 2& ja

156

p,(r)

22 August 1986

CHEMICAL PHYSICS LETTERS

Volume 129, number 2

t h,

(2)

where a = a0 t 612, the distance of closest approach of the small ions to the polyion axis; and pc(r) is the counterion radial distribution function. After the MC simulation has been equilibrated, /3(A) is evaluated numerically by summing the histogram representation of p,_(r) over the appropriate number of radial increments. More details concerning the calculation of pe(r) and /3(A) are given in the appendix of paper I.

3. Results and discussion For counterions of diameter 6 = 0.6 nm, MC radial distributions were computed for both the helical and the axial model of the DNA charge distribution at several salt concentrations in the range 1.5 < pa]/ [p] < 7.0. These distributions were used in eq. (2) to calculate the values of P(A) that are shown in fig. 1. Also shown are the corresponding predictions of the PB equation and pee , the number of condensed ions predicted by CC theory. Fig. 1 indicates that, for an ion diameter of 0.6 mn, flhax(A) > /3”“(A) > flpB(A). The increase of /3line(A) over pm(A) results from the proper treatment of small ion correlations in the MC simulations, as discussed previously [3,4,6,7]. The dis-

2.0

4.0

6.0

8.0

No/P Fig. 1. Values of o(O.3) and p(O.7) as functions of salt concentration ([Na]/[P] with [P] = 0.031 mol/dmj) for an ion diameter of 0.6 nm. The MC value8 of @are evaluated with eq. (2) from simulations using either a helical model for DNA (o), or a linear model (0). The PB values (a) were calculated using numerical solutions of the cylindrical PB equation. The salt-invariant value of flee (the number of condensed ions per DNA phosphate predicted by the CC theory) is also shown I-).

Volume 129, number 2

CHEMICAL PHYSICS LETTERS

Crete helical model charge distribution produces a 5--10% increase in /3hen”(A) over fin”(A). Simulations of small-ion distributions in the vicinity of a plane having a uniform discrete surface charge density also show enhancements of p&‘(A) over /Icon(A), computed near a charged plane with a comparable continuous surface charge density [8,9]. According to Manning’s molecular CC theory for DNA in a uniunivalent salt solution [ 1 I], 0.76 counterions per phosphate charge are located within the annular volume, VP, extending from the unhydrated radius of DNA, ao = 1.Orun, to ml .7 nm. To compare flee = 0.76 with the MC values of p(A), it is necessary to select the volume of thickness A that corresponds to the VP parameter, which is evaluated in CC theory without considering the finite size of small ions. Since realistic ion sizes are assumed in the MC simulations, the MC annular volume corresponding to VP is chosen to extend outward 0.7 nm from the distance of closest approach of the small ions to the DNA cylinder (1.3 run). Thus, p(O.7) is the number of counterions per polyion charge whose centers are contained in the adjacent annular volume that is the same size as VI,. Over the entire range of salt concentrations inveatigated here, 0” is considerably greater than @h*x(0.3), but @hefix(0.7)brackets the CC value (phGx(0.7) = 0.58 at [Na]/P] = 1.5 and phenx(0.7) = 0.78 at pal/P] = 7.0). For both A = 0.3 nm and A = 0.7 nm, /IhelIx displays essentially the same salt dependence as do pnne(A) and flPB(A). In contrast, the CC theory predicts that the number of condensed ions is (nearly) salt-invariant under these conditions. For the conditions and model assumptions indicated above, the presence of a discrete helical array of charges on the surface of the model polyion does not significantly alter the salt dependence of P(A) that is predicted by MC and PB calculations for DNA modeled with an axial charge density [4]. The MC results discussed thus far were obtained for small ions of diameter 6 = 0.6 nm. Since ion distributions calculated with the PB equation are known to be sensitive to the choice of a, it is of interest to test whether p(A) is affected by decreasing the distance of closest approach of the mobile ions to the polyion charges. For this purpose some simulations were conducted for small ions of diameter 6 = 0.1 nm. Fig. 2 demonstrates that a reduction in the ionic diameter to 0.1 run in the MC simulations does not alter the rela-

22 August 1986

0.5 0.4

Fig. 2. Values of p(0.3) and p(0.7) as functions of salt concentration for an ion diameter of 0.1 mu. The symbols de noting the values of p calculated in different ways are the same as in fv. 1.

of @hre(A) over pPB(A) calculated for the case of 6 = 0.6 mn. However, for 6 = 0.1 nm . the values of /Ihelix change radically in comparison to P”““(A), as shown in fig. 2 for flhenx(0.3) and /Ihenx(0.7) at three salt concentrations. A substantial enhancement in /Ihe1ix(A)is observed over pEne(A) and pPB(A). The magnitude of /?e1ix(0.3) is closer to PC’ than is either /Iline(0.3) or flw(0.3). This striking increase in phenx(A) results from formation of contact ion pairs between the smaller counterions and the discrete polyion charges (observed in the MC histograms as the immobilization of counterions in the near vicinity of polyion charges). Thus, the pairing of small ions (8 5 0.1 nm) with polyion charges can produce a substantial increase in the magnitude of flheEx(A) over the value calculated for 6 = 0.6 run; however, this ion pairing does not significantly alter the salt dependence of /Ihelix( The ion pairing observed in our MC simulations for (unrealistically) small ions potentially introduces an additional uncertainty into MC estimates of /Ihe’“( In this situation Monte Carlo calculations with all ions initially free in solution are exceedingly time consuming because of the extremely slow convergence of /3henx(0.1). Occasionally, a small counterion randomly samples the region in the near vicinity of a polyion charge and becomes essentially immobilized. This “ion pairing” then produces a sudden increase in /3hehx(0.1). To test whether /Ihelrx(0.1) had converged under these tive enhancement

157

Volume 129, number 2

CHEMICAL PHYSICS LETTERS

circumstances, the effect of variations in me mitial specification of small ion coordinates was investigated. Simulations in which some counterions initially were paired with polyion charges were compared with simulations in which all ions were initially free in solution. Both sets of simulations were used in computing the average values of /Ihehx(A) reported in fig. 2. Typically, when a simulation is initiated with no ions paired, the surface concentration, ~hehx(O.l)/rrb [(a + 0.1)2 - a2], increases from ~3.0 mol/dm3 very slowly and finally converges at a value in the range 8-9 mol/dm3 after 2 X lo6 configurations. However, when a simulation is initiated with several ions paired (so that the initial surface concentration is in excess of 10 mol/ dm3) the surface concentration falls to 8-9 mol/dm3 in (S-10) X lo5 configurations. Consequently, the choice of the initial set of coordinates can dramatically affect the convergence rate of /3he1ix(0.1),but it ap parently does not significantly alter the value at which convergence is achieved.

4. Conclusions For a counterion diameter S = 0.6 run, use of a discrete helical charge model in place of the axial charge model of DNA yields values of flh*“(A) that are somewhat greater than @““(A), but does not affect the salt dependence of /3(A) for any value of A over the range of salt concentrations investigated (1.5 < ma]/ [p] < 7.0). For a physically unrealistic counterion diameter (8 = 0.1 mn), counterion immobilization occurs at some polyion charges during the MC simulations, and thus p(A) substantially increases. Although the resultant magnitudes of p(A) are in better agreement with /ICC, the salt dependence of flhe”“(A) contrasts both with the prediction of CC theory and with inferences from NMR data. (These NMR studies also provide no experimental evidence for contact ion pairing.) The substantial salt dependences exhibited by both /Ihelix and /Ifine indicate that more details (possibly including counterion binding, local dielectric effects, the local structure and penetrability of the DNA helix, for example) must be built into the model used in the MC computations in order to account for the salt invariance of the number of counterions close to DNA (and other highly charged

158

22 August 1986

polymers) that has been inferred from NMR measurements. It is also of interest to understand why the molecular and thermodynamic predictions of CC theory appear to be in better agreement with experiment than the simplicity of its model should allow.

Acknowledgement This research was supported by NIH Grant GM3435 1. Computations were performed at the IBM Palo Alto Scientific Center (PASC) Computing System through a generous grant of computing time to the Scientific Research Program at the Department of Chemistry, UW Madison. We are grateful to A. Karp and R.A. Blaine of the PASC for their assistance in implementing this grant.

References [l] D. Bratko and V. Vlachy, Chem. Phys. Letters 90 (1982) 434. [2] M. LeBret and B.H. Zimm, Biopolymers 23 (1984) 271. [ 31 C.S. Murthy, R.J. Bacquet and P.J. Rossky, J. Phys. Chem. 89 (1985) 701. [4] P.A. Mills, C.F. Anderson and M.T. Record Jr., J. Phys. Chem. 89 (1985) 2984. [S] P. Lime and B. Jonsson, J. Chem. Phys. 78 (1983) 3167. [6] B. Jiinsson, H. Wennerstriim and B. Halle, J. Phys. Chem. 84 (1980) 2179. [7] GM. Torrie and J.P. Valleau, J. Chem. Phys. 73 (1980) 5807. [8] W. van Megen and I. Snook, J. Chem. Phys. 73 (1980) 4656. [9] I. Snook and W. van Megen, J. Chem. Phys. 75 (1981) 4104. [ 101 G.S. Manning, Quart. Rev. Biophys. 11 (1978) 179. [ll] G.S. Manning, Accounts Chem. Res. 12 (1979) 443. [ 121 J. Reuben, M. Shporer and E. Gabbay, Proc. Natl. Acad. Sci. US 72 (1975) 245. [13] L. Herwats, P. Laszlo and P. Gerard, Nouv. J. Chhn. 1 (1977) 173. [14] C.F. Anderson, M.T. Record Jr. and P.A. Hart, Biophys. Chem. 7 (1978) 301. [ 151 M.L. Bleam, C.F. Anderson and M.T. Record Jr., Biochemistry 22 (1983) 5418. [16] L. Nordenskiold, D.K. Chang, C.F. Anderson and M.T. Record Jr., Biochemistry 23 (1984) 4309. [ 171 W.H. Braunhn and L. Nordenskiold, Eur. J. Biochem. 142 (1984) 133. [18] B.J. Klein and G.R. Pack, Biopolymers 22 (1983) 2331.