J. Phys. Chem. B 2008, 112, 13359–13366
13359
Estimation of the Persistence Length of DNA from the Torsion Elastic Constant and Supercoiling Free Energy: Effect of Ethylene Glycol David P. Rangel, Bryant S. Fujimoto, and J. Michael Schurr* Department of Chemistry, UniVersity of Washington, P.O. Box 351700, Seattle, Washington 98195-1700 ReceiVed: February 20, 2008; ReVised Manuscript ReceiVed: April 25, 2008
Two different methods are proposed to estimate the persistence length (P) of DNA from the measured torsion elastic constant (R) and the twist energy parameter (ET) that governs the supercoiling free energy. The first method involves Monte Carlo simulations and reversible-work calculations of ET for model DNAs that possess the measured R and selected trial values of P. Comparison of the computed ET values with the experimental value allows estimation of P (or equivalently the bending elastic constant (κβ)) by interpolation. A far simpler, though less accurate, alternative is to solve a previously conjectured analytical relation connecting ET, R, κβ (or P), and an unknown “constant” (B). The present simulations are used to ascertain the optimum value of B and to assess the validity and accuracy of that relation. Within the simulation errors, P values obtained from the measured R and ET via this analytical expression agree with those determined from the simulations and ET values reckoned from the input R and κβ by this analytical expression agree with the corresponding simulated values. Although B is found to be insensitive to variation in R, it appears to decline slightly with increasing κβ. The original analytical expression is modified to take this apparent variation of B with κβ into account. By using this modified analytical relation to estimate P (from the measured R and ET) or ET (from the input R and κβ), much closer agreement is obtained respectively with the values of P or ET obtained from the simulations. As specific examples, these methods are applied to determine P in 0 and 20 w/v % ethylene glycol, which has been shown to induce a structural transition in duplex DNA. Introduction Numerous conformational substates of duplex DNAs within the right-handed B and A families have been detected by various methods.1 Properties that are sensitive to changes in secondary structure respond to practically every conceivable perturbation, including the addition of small neutral osmolytes.1-4 Such conformation-sensitive properties include the elastic constants, R and κβ, for torsion and bending, respectively, and the twist energy parameter, ET, that governs the supercoiling free energies of circular topoisomers. Measurements of R by time-resolved fluorescence polarization anisotropy (FPA) and ET by the equilibrium topoisomer distribution method can be performed under conditions, where any imposed deformational strain is either nonexistent or so slight that it does not significantly alter the equilibrium between conformational substates. Consequently, any changes in R and κβ arise exclusively from changes in intensive variables (e.g., T and concentrations of various ligands, neutral osmolytes, or counterion species). It is also desired to assess any changes in κβ, which could in principle be obtained by a measurement of the persistence length using some method that does not involve substantial imposed bends or tensile stresses. As will be shown, κβ can also be determined from the measured R and ET by computationally intensive simulations or estimated more approximately from R and ET by using either of two simple analytical expressions containing unknown constants that are evaluated in the present study. Effects of the small neutral osmolyte, ethylene glycol (EG), on the supercoiling free energy and intrinsic twist of p30δ DNA were recently investigated.2-5 Increasing the EG concentration from 0 to 40 w/v % at 310 K caused a decreasing (or reverse) * Corresponding author. E-mail:
[email protected]. Telephone: 206-543-6681. Fax: 206-685-8665.
sigmoidal transition in the graph of ET versus -ln aw, where ET is the twist energy parameter that governs the supercoiling free energy and aw is the activity of water. EG had a similar effect on the intrinsic twist. The midpoint of this transition occurs near 20 w/v % EG and its final plateau is reached at ∼30 w/v % EG. The full profile of this transition will be presented in a subsequent paper. Increasing the EG concentration from 0 to 19.4 w/v % caused a 1.15-fold decrease in ET, a 1.27- to 1.35-fold increase in the torsion elastic constant between base pairs, R, and a modest but statistically significant change in the circular dichroism (CD) spectrum.4,5 The bending elastic constant was not directly measured. Nevertheless, from the observed decrease in ET and increase in R, one may infer qualitatively that the average equilibrium bending elastic constant between base pairs, κβ, must undergo a significant decrease. In the present work, quantitative estimates of the κβ values in 0 and 19.4 w/v % EG are obtained from the corresponding measured values of ET and R by Monte Carlo (MC) simulations and reversible work calculations of the supercoiling free energy. This protocol provides a novel way to estimate κβ, which in turn yields the corresponding values of the equilibrium bending rigidity, A ≡ hκβ, where h is the rise per base pair, and yields the equilibrium persistence length, P ) hκβ/kT, where k is Boltzmann’s constant and T is the absolute temperature. The miniscule EG-induced change in intrinsic twist (∼0.4 turn out of 474 turns) and modest change in the CD spectrum suggest that the DNA remains in the B family and that changes in the structure are rather small. For simplicity, 19.4 w/v % EG will henceforth be referred to simply as 20 w/v % EG. A simple analytical expression connecting ET to R and κβ would enable one to evaluate any one of these three constants
10.1021/jp801526w CCC: $40.75 2008 American Chemical Society Published on Web 08/22/2008
13360 J. Phys. Chem. B, Vol. 112, No. 42, 2008 from independent estimates of the other two, without resorting to computer simulations in every case. Such a relation was derived previously by assuming that (1) the distribution of writhe in a long nicked circular DNA is Gaussian, which implies a quadratic potential of mean force for writhing, and (2) the Hooke’s law torque constant for writhing, κw, is proportional to κβ with proportionality constant B, so κw) Bκβ.6 Neither the validity of the resulting expression for ET (eq 3 below) nor the value of B has been investigated previously for any model of a long DNA with a realistic repulsive potential. A second objective of the present study is to test eq 3 for ET and to ascertain the value of B under standard conditions by means of MC simulations. We also examine the sensitivity of B to the change in torsion elastic constant, R, that occurs upon changing from 0 to 20 w/v % EG, as well as its possible variation with changes in bending elastic constant, κβ. Values of the average torsion elastic constant, R, were measured for linear p30δ DNA under various conditions via the time-resolved fluorescence polarization anisotropy of intercalated ethidium (i.e., the FPA method).2-7 The measured value in 0 w/v % EG, namely, R ) 4.99 × 10-12 dyne cm, lies in the normal range for other linear and weakly strained DNAs under the same conditions. Substantially larger R values have been obtained for small (181 e N e 250 bp) circular DNAs, which are highly bent (g1.5 °/bp), by topoisomer ratio (TR), cyclization kinetics (CK), and FPA methods, and for much longer DNAs under considerable tension by single-molecule methods, as noted previously.8 However, recent MC simulations of the supercoiling free energy indicate that such large R values do not apply to weakly strained, large (N g 2600), circular DNAs with small linking differences.8 In contrast, MC simulations with the canonical persistence length, P ) 500 Å, and R values obtained for linear and weakly strained large circular DNAs by the FPA method yielded supercoiling free energies in good agreement with experiment.8 Moreover, when the FPA method was applied first to a linear and then to a circularized form of the same 181 bp DNA, the torsion elastic constant jumped from a value typical of linear and weakly strained DNAs to higher values typical of TR and CK measurements on small (N e 215 bp) circular DNAs.9 In addition, when the bending strain was reduced by increasing the size of small circles to 350-360 bp, the R value obtained by the CK method fell into the range observed for linear and weakly strained DNAs.8 These and other observations imply that the FPA method yields the correct R value for a given DNA with its particular mean deformational strain and suggest that sufficient bending strain or tensile stress shifts the conformational equilibrium in favor of substates with a greater average torsion elastic constant relative to that of relaxed linear and weakly strained DNAs. In summary, R values obtained by the FPA method are appropriate for the present simulations of circular p30δ DNAs with small linking differences. The simulation protocol adopted here is similar to that described previously,10,11 but the intersubunit repulsive potential is that proposed by Delrow et al.,12 which provides a more realistic account of the electrostatic interactions. When the value of R measured for p30δ DNA and the canonical persistence length, P ) 500 Å, were used as input parameters, the previously simulated values of ET for p30δ were in good agreement with the measured values over the entire range of superhelix density from 0 to native,10 and the predicted shapes of equilibrated surface-confined plasmids were in good agreement with AFM images.11,13 A more general protocol was devised to simulate model DNAs with nonuniform elastic constants and, when applied to partially melted supercoiled DNAs, yielded remark-
Rangel et al. ably good agreement with certain experimental results pertaining to ET.14,15 That same protocol also yielded good agreement with the ET values measured for p30δ and many other DNAs at low superhelix density, provided that a torsion elastic constant typical of unstrained or only weakly strained DNAs was employed in the simulations.8 Plan of the Paper Our first objective is to simulate a model p30δ DNA with small linking difference in 0 w/v % EG using the measured value of R and different trial values of the persistence length, P, and then to compute ET in each case by a reversible work protocol. The “correct” value of P is determined from a comparison of the computed ET values with the experimental ET datum via an interpolation protocol devised for that purpose. This entire procedure is reiterated for p30δ DNA in 20 w/v % EG, and the “correct” value of P in that solution is similarly determined. Our second main objective is to test the validity of a previously proposed simple analytical expression for ET in terms of R and κβ, to evaluate an unknown “constant” (B) therein and to assess the variation (if any) of B with R and κβ. The following sections present (i) the basic theory, (ii) the mesoscopic model, (iii) the simulation protocol and input parameters, (iv) the calculation of the supercoiling free energy, (v) certain details of the computations, (vi) the presentation and discussion of the results, and (vii) the conclusions. The interpolation protocol is described in Results and Discussion. A subsequent section describes (i) testing of the analytical formula relating ET to R and κβ, (ii) evaluation of the “constant” B therein, and (iii) assessment of the possible variation of B with R and κβ. Theory For a given topoisomer, the number of turns, l, of one strand around the other is a topological invariant called the linking number. The state of deformation of a supercoiled DNA is characterized by its linking difference, ∆l ) l - l0, where l0 ) Nφ0 is the nonintegral intrinsic twist and φ0 is the intrinsic succession angle between base pairs. The linking number is partitioned between the twist (t) and the writhe (w) according to l ) t + w,16,17 so
∆l ) t - l0 + w
(1)
Both simulations8,10,14 and experiments2-4,18 on long DNAs (N g 2600 bp) indicate that the change in supercoiling free energy to alter the linking difference from ∆m to ∆l in ∼0.1 M NaCl is accurately represented by
∆Gsc ) kT(ET /N)(∆l2 - ∆m2)
(2)
wherein ET is independent of ∆l and ∆m, especially over the thermally populated range of topoisomers. The previously derived analytical expression for ET in terms of the elastic constants, R and κβ, is given by,6
ET ) ((2π)2/(2kT))(RBκβ/(R + Bκβ))
(3)
where B is the unknown proportionality constant in the assumed relation, κw ) Bκβ. After reckoning ET for different values of κβ and R under the same solution conditions via Monte Carlo simulations, we can solve eq 3 for the value of B in each case and determine its variation, if any, with variations in κβ and R.
Estimation of the Persistence Length of DNA
J. Phys. Chem. B, Vol. 112, No. 42, 2008 13361
Mesoscopic Model The circular DNA containing N base pairs (bp) is replaced by a circular chain of M rigid rod subunits, each containing ν ) N/M bp. These subunit rods are labeled consecutively by the index j, j ) 1,..., M. In each subunit (j) is fixed a coordinate frame (xj,yj,zj), the zj axis of which lies along the bond vector (bj) that extends from the origin of the jth frame to the origin of the succeeding (j + 1)th frame. The Euler rotation that carries a coordinate frame from coincidence with the jth frame to coincidence with the (j + 1)th frame is Φj,j+1 ) (Rj,j+1βj,j+1γj,j+1), where the component rotations, Rj,j+1, βj,j+1, and γj,j+1 are taken sequentially around the body-fixed zj, new body fixed y′j, and final z′′j axes.19 The net twist (turns) of the DNA is given by M
t)
∑ φj,j+1/2π
(4)
j)1
where φj,j+1 ≡ Rj,j+1 + γj,j+1 (rad) is the net twist of the Euler rotation from the jth to (j + 1)th frame.6,8,10,11,14,20 The writhe is approximated by the discretized Gauss integral21
w)
(
1 4π
M
M
)∑ ∑
i)1 j)1j*i
(bj × eij · bj) |ri - rj|2
(5)
wherein ri and rj denote the positions of the origins of the ith and jth subunit frames, respectively, in the laboratory frame, and eij ) (ri - rj)/|ri - rj| denotes a unit vector along ri - rj. Each rigid rod subunit of the model DNA is connected to its neighbors at either end by Hookean torsion and bending springs. The subunit length, b ) 95.4 Å, was chosen to be less than 1/5 of a persistence length, P ∼ 500 Å. Various properties, including the mean squared writhe of a nicked circular DNA and the supercoiling free energy of a closed circular DNA, were found to become independent of rod length, when b e (P/5), as shown by simulations of nicked circles22 and also by (unpublished) simulations of closed circles in our laboratory. This choice of b corresponds to ν ) 28.06 bp per subunit. The number of subunits was chosen to be M ) 175, which corresponds to N ) 4910 bp, close to the actual number, N ) 4932 bp, determined by direct sequencing of p30δ.4 The input values of the torsion elastic constants are taken to be the measured values, corrected for the slightly lower temperature and higher Mg2+ concentration prevailing in the measurements of ET, as described subsequently. Simulations were performed using different trial values of the persistence length, as detailed below. Comparison of the computed ET values with the experimental data then allows interpolation of the prevailing value of P in each case. Potential Functions The potential of mean force of any particular configuration is given by
Utot ) Utw + Ube + UI
(6)
where Utw and Ube are the torsion and bending potential energies, respectively, and UI ) Uhc + Uel is the intersubunit potential energy, which includes both hard-cylinder (Uhc) and electrostatic interactions (Uel) between different subunits of the model DNA.8,11,12,14 Measurements of the kinetics of cyclization of 200 bp circles with a designed intrinsically straight sequence and also circles with a native sequence suggest that any sequence-dependent intrinsic bends in the native sequence make a negligible contribution to the curvature in small circles.23 It is assumed here, as in most of our previous work, that sequence-dependent intrinsic bends
contribute negligibly to the mean-squared curvature and persistence length of native DNA. The potential energy function for torsion is given by6,10
Utw )
M ˜ ˜ R R(2π)2 (φj,j+1 - φ0)2 ) (∆l - w)2 + 2 j)1 2M
∑
˜M R (φ - 2πt ⁄ M)2 (7) 2 j)1 j,j+1
∑
where R˜ ≡ R/ν ) Rh/b is the elastic constant of the torsion springs between rigid subunits, which is chosen to yield the same overall torsional rigidity as that of a filament with torque constant R for the springs between base pairs. In the last line, the first term is the contribution of the mean net twisting strain per subunit, (2π)(t - l0)/M (rad) and represents the minimum energy required to introduce the net twist t into a molecule with a fixed configuration of its bond vectors. The second term is due to fluctuations of the local twist about that mean net value. Because UI is assumed to be independent of the azimuthal orientations of the interacting subunit rods, it is independent of such fluctuations in local twist. Consequently, the contribution of these local twist fluctuations to the partition function can be evaluated exactly by analytical methods and is always the same quantity, independent of t, w, l, or knot topology.6,10 Hence, this second term in eq 7 makes no contribution to the relative probabilities of different configurations of the bond vectors or to the difference in free energy between any two topoisomers. Consequently, for the purposes of reckoning supercoiling free energies or ensemble-averaged configurational properties, the second term in eq 7 can be omitted. However, if the elastic constants are not uniform and isotropic, a more general protocol must be employed.8,14,20 The subunit rods are regarded as hard cylinders of 24 Å diameter. Uhc is infinite if the cylinders overlap, and it is zero otherwise. The algorithm to check for overlap was described previously.11 The electrostatic potential energy, Uel, is modeled as follows.8,11,12,14 Three charged spheres of radius a ) 12 Å are spaced every 31.8 Å along each bond vector, beginning with one at the tail-end of each, which coincides with the head-end of the previous bond vector. The electrostatic potential of mean force between the ith and jth charged spheres is given in cgs units by12
Uel )
(Ze)2 exp[-λa] ε 1 + λa
(
)( 2
exp[-λrij] rij
)
(8)
where λ is the Debye screening parameter, e is the protonic charge (esu), ε is the relative dielectric constant, and rij is the distance between the ith and jth charged spheres. The effective valence of the sphere is -Z and is chosen in the following way. First, the nonlinear Poisson-Boltzmann (NLPB) equation is solved for an infinite cylinder of 24 Å radius with the linear charge density of DNA in the appropriate ionic medium. Each charged sphere is assumed to exhibit a Debye-Hu¨ckel screened Coulomb electrostatic potential (esu) with effective charge, -Ze. The sum of Debye-Hu¨ckel potentials of 2001 spheres in a straight DNA is evaluated along each of a series of lines perpendicular to the chain near its center, and the average potential at a given distance, r, from the chain of spheres is evaluated. Z is adjusted so that this average potential closely matches the NLPB result at all distances, r g 24 Å. Excellent matches were achieved in the present work, as found also in previous studies.8,11,12,14 The electrostatic interaction energy in
13362 J. Phys. Chem. B, Vol. 112, No. 42, 2008
Rangel et al.
TABLE 1: Simulated Total Persistence Length, Ptot, of Linear Chains Comprising 10, 20, and 30 Subunit Rods of Length 95.4 Å in 0 and 20 w/v % Ethylene Glycol (EG) w/v % EG
Pbe
N
〈Ptot〉
σ〈Ptot〉
〈〈Ptot〉〉a
0 0 0 0 0 0 20 20 20 20 20 20 20 20 20
300 300 300 500 500 500 284 284 284 300 300 300 500 500 500
10 20 30 10 20 30 10 20 30 10 20 30 10 20 30
304 298 300 505 504 497 274 284 287 290 293 301 485 494 500
12 8 6 21 15 12 10 8 6 13 7 7 20 15 13
300 ( 3
Z
λD (Å)
ε
R × 1012(erg)
P (Å)
0 20
-10.91 -10.54
8.7790 8.4472
74.0 68.2
6.308 8.536
300, 500 284, 300, 500
Simulation Protocol 282 ( 7 296 ( 5 493 ( 7
eq 8 is the standard Derjaguin-Landau-Verwey-Overbeek expression for the interaction between two charged spheres. In the present study, the electrostatic interaction energy is evaluated between all pairs of spheres that are separated by two or more adjacent spheres along the chain but which also are less that 220 Å apart. The former requirement practically eliminates the contribution of electrostatic repulsions between adjacent subunits to the bending rigidity and persistence length. The latter requirement neglects all electrostatic interactions that are
