INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS GLOSSARY REFERENCES

The large-scale fluctuations of DNA are key to understanding kinetically complicated events, such as the ease of the long, threadlike molecule snaking though the pores of a gel or closing into a loop between separately bound regulatory proteins. DNA loop formation is implicated, in turn, in a number of important biological processes, including the regulation of transcription (Bellomy et al., 1988; Schleif, 1992) and the organization of chromatin (Dillon et al., 1997; Bazett-Jones et al., 1999; Ringrose et al., 1999).

The motions of polymeric DNA are generally modeled in terms of a spatially homogeneous, naturally straight, elastic rod which ignores realistic features of chemical structure (reviewed in Olson, 1996). The fluctuations and correlations of structural parameters in crystals of pure DNA and of DNA-protein complexes, however, show that the equilibrium rest states and elastic constants of neighboring base pairs are dependent on sequence and further reveal the presence of strong couplings between modes of deformatione.g., between bending, twisting, and stretchingthat are the direct result of the chiral nature of DNA chemical structure (Olson et al., 1998). These local features become important as the deformations of sequential base-pair steps accumulate with the increase of chain length, introducing notable structural variability at the mesoscopic level, i.e., in chains of a few hundred base pairs.

Electron micrographs and other low resolution images of mesoscopic DNA fragments reveal broad distributions of molecular shapes and end-to-end distances (Muzard et al., 1990; Thresher and Griffith, 1992; Bednar et al., 1995; Lyubchenko et al., 1995; Hansma et al., 1996; Bustamante et al., 1997). The sequential features responsible for these spatial arrangements are often deduced by comparison of the experimentally observed images of a given DNA with the distribution of configurations generated by Metropolis-Monte Carlo sampling of the likely structural fluctuations of the constituent dimers (Hagerman, 1985; Levene and Crothers, 1986; Kahn and Crothers, 1998). Limitations on computer resources make routine comparison of simulations and experiment impractical, because hundreds of thousands of coordinate sets must be generated to characterize an individual polymer molecule. The study of the large-scale, collective motions of long DNA is consequently impeded.

The normal modes of DNA provide an alternative, computationally less demanding way to study the equilibrium properties of DNA. The normal modes are coupled vibrations found by assuming that the molecular potential energy can be approximated by a harmonic function of the configurational variables and by then solving a generalized eigenvector problem to give a closed analytical description of the motion. The eigenvalues give the vibrational time scales (frequencies) and the eigenvectors the details of the corresponding motion. The motion can be described and visualized at each frequency or as a superposition of independent modes. Time-averaged equilibrium and kinetic properties can be calculated accurately and efficiently as weighted sums, e.g., thermodynamic functions can be determined from the normal mode frequencies. The price paid for these benefits is the limited accuracy of the harmonic energy approximation. The method, nevertheless, provides a useful first impression of the flexibility of a molecule and shows how the motions may change when the chemical structure is modified or the chain is perturbed by the binding of proteins and other ligands.

Normal mode or vibrational analysis is a well established and much used technique for the study of both small molecules and proteins: the technique is successful in reproducing the vibrational spectra of small molecules (Wilson et al., 1955) and in analyzing the collective motions involved in the folding of protein fragments and in the binding and release of substrates and products to and from enzymes (Levy and Karplus, 1979; Brooks and Karplus, 1983; Kitao and G, 1999; Berendsen and Hayward, 2000). Analyses of DNA normal modes to date have focused on the collective motions of very short chain fragments ranging from small oligomers (Tidor et al., 1983; Irikura et al., 1985; Garcia and Soumpasis, 1989; Kottalam and Case, 1990; Ha Duong and Zakrzewska, 1997a,b, 1998; Lin et al., 1997) to a few turns of double helix (Matsumoto and G, 1999). The size of the latter systems is limited by the kinds of parameters which have been used to describe three-dimensional structuretypically the Cartesian coordinates of the constituent atoms or the dihedral angles of a nucleic acid fragment with rigid chemical bonds and valence angles. The small, almost imperceptible motions of short DNA oligonucleotides described at this high level of resolution are closely tied to the detailed sequence-dependent fine structure of the double helix and the association of water, drugs, proteins, and other molecules with individual nucleotide residues. The conformational effects which are collectively responsible for the large-scale, sequence-dependent properties of polymeric DNA chains are difficult to discern in such studies. The normal modes of infinitely long DNA, however, have been examined by Prohofsky and associates (Hua and Prohofsky, 1988; Chen and Prohofsky, 1995) by repeating a short fragment of DNA and taking advantage of helical symmetry to reduce the number of independent coordinates. The complexity of chains that can be studied from this perspective is again limited by the length of the repeating unit.

In the present investigation, we examine the motions of long fragments of DNA by taking advantage of "knowledge-based" harmonic energy functions deduced from known high-resolution crystal structures of DNA (Olson et al., 1998). The deformations of individual dimers are described by six independent "step" parameters which specify the spatial arrangements of neighboring base pairs: three angular variables called Tilt, Roll, and Twist and three variables called Shift, Slide, and Rise with dimensions of distance (Dickerson et al., 1989). This rigid-body representation of base pairs significantly reduces the number of independent variables per chain molecule, thereby making it possible to study the normal modes of much longer DNA fragments than could heretofore be treated in either Cartesian or dihedral angle space. In addition, the elastic energy incorporates the sequence-dependent anisotropy of DNA bending as well as the known correlations of base-pair step parameters. The only missing structural information is the detailed conformation of the sugar-phosphate backbone, including the charged phosphate groups. The latter atoms and the surrounding aqueous solvent, such as the water molecules and counterions in the first or second solvation layer around the DNA, are implicitly treated in the energy terms so that their omission introduces no serious errors when duplex deformations are limited to energies of the order of k_BT (where k_B is the Boltzmann constant and T the temperature in Kelvin). It should be noted, however, that the surrounding solvent does not have a viscosity and that we do not consider the damping effect of solvent on large-scale polymeric properties in this study. This low-resolution model, nevertheless, provides a straightforward way to deduce the effects of sequence on global folding, something that is said (Berendsen and Hayward, 2000) still to elude analyses of the slow conformational changes of proteins. In principle, there is no limit on the length of DNA that can be described in this manner so long as the constituent dimers obey the harmonic energy model, i.e., there are no excursions of the molecule from the classical B-form double helical structure. Longer chains are characterized by a greater number of normal modes, which when appropriately weighted and superimposed, give rise to a broader distribution of global molecular configurations. In practice, chain length is restricted by the form of the kinetic energy, which assumes that changes in atomic coordinates brought about by the fluctuations of individual "step" parameters are small. These upper limits preclude the need to treat the long-range electrostatic self-repulsion that influences the folding of long supercoiled DNA molecules (Fenley et al., 1994; Vologodskii and Cozzarelli, 1995; Westcott et al., 1997).

We focus attention here on the normal modes of regularly repeating, linear polymers free of bound ligands to establish a point of reference for studies to be reported elsewhere of arbitrary DNA sequences and of chains with ends held in place by specific anchoring conditions, such as the looped configurations imposed by protein binding. We compare our simplified energy model with previous all-atom treatments of short oligomers and with ideal elastic rod representations of polymeric DNA. We identify the local base-pair fluctuations responsible for the dominant (lowest frequency) normal modes. The derived polymeric properties account satisfactorily for the overall bending, i.e., persistence length, of the double helix and reveal the critical role of sequence in global twisting. The computations further suggest that the constants impeding the large-scale stretching of single DNA molecules are related to the concerted bending, twisting, and lateral displacements of neighboring base-pair planes as well as to their axial separation. We identify particular sequence contexts under which ideal elastic behavior breaks down, concentrating on conditions that induce the mesoscopic bending anisotropy which underlies DNA loop formation. Finally, by comparing the end-to-end dimensions of the simulated double helices with exact values obtained from standard matrix formulations of polymer configurational statistics (Flory, 1969; Maroun and Olson, 1988; Marky and Olson, 1994), we investigate the chain lengths at which the normal mode analysis of DNA modeled in terms of base-pair steps is valid.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS GLOSSARY REFERENCES

Base-pair representation

The atoms of each base pair lie in the plane of a local, orthogonal coordinate frame (x, y, z), which is defined in accordance with recently established guidelines (Olson et al., 2001). (Note: The x and y axes lie in the plane of the base pair, with x pointing in the direction of the major groove along the pseudodyad axis of the base pair and y running along the long axis of the base pair in the direction of the leading (sequence) strand, parallel to the C1' ··· C1' vector, and displaced so as to pass through the intersection of x with the vector connecting the pyrimidine Y(C6) and purine R(C8) atoms. The z axis is perpendicular to the base-pair plane, pointing in the 5'-3' direction of the leading strand.) The coordinate system of the kth base pair is further characterized by the position of its origin o_k and the rotation matrix R_k, which relates the local base-pair frame to a fixed, orthonormal global reference frame. The rotation matrix is given by [_xyz], where the ( = x, y, z) are unit vectors, expressed in columns, along the positive coordinate axes of the base-pair frame.

The relative position and orientation of the (k + 1)th base pair with respect to its predecessor k are given respectively by the difference between coordinate origins, v_k,k+1 = (o_k+1  o_k) and the product of the rotation matrix R_k+1 and the inverse of R_k, i.e., the dimer transformation matrix T_k,k+1 = R_k+1(R_k)¹. The components of the projection of v_k,k+1 on the coordinate axes of base pair k, (_x, _y, _z), are used as translational parameters in the present calculation.

The matrix T can also be expressed in terms of the rotation of magnitude  around the axis u which brings the coordinate frames on successive base pairs into coincidence, where the components of the unit vector u are (u_x, u_y, u_z) and the angle  equals ( +  + )^1/2 (Chasles, 1830). The elements of T, so expressed, are given by (Jeffreys and Jeffreys, 1946):

(1)

where _µ is the Kronecker delta, i.e., _µ = 1 when  = µ, _µ = 0 when   µ, and _µ = ±1 when , µ,  is an even or an odd permutation of 1, 2, 3, respectively, and vanishes otherwise. The rotational components (_x, _y, _z) are equated to (Tilt, Roll, Twist) after the definition of Babcock et al. (1994).

The translational parameters between base-pair planes are generally expressed in terms of a "middle" frame so that numerical values are independent of the direction from which a DNA structure is analyzed (Lu and Olson, 1998), i.e., either the sequence or the complementary strand. The translational parameters defined in this work in terms of the kth coordinate frame are related to Shift, Slide, Rise values in the literature, e.g., (Olson et al., 1998), through the matrix T^1/2 that effects the halfway rotation (/2) between neighboring base-pair planes:

(2)

The dependence of the ( = x, y, z) on Tilt, Roll, Twist (through their incorporation in T) as well as on Shift, Slide, Rise precludes description of the local conformational energy in terms of a quadratic expression in the six base-pair step parameters. To overcome this limitation, we express the rotation matrix in Eq. 2 in terms of the equilibrium "rest" values of Tilt, Roll, and Twist, i.e., T^1/2  T_o^1/2. The latter approximation is valid so long as the dimer structure does not depart significantly from this minimum energy state.

Normal mode analysis

Our treatment of the normal modes of DNA at the level of "step" parameters builds upon general formulations (Noguti and G, 1983; Braun et al., 1984; Higo et al., 1985; Levitt et al., 1985) previously developed to describe the collective motions of proteins in terms of internal chemical coordinates, i.e., dihedral angles. The conformational potential energy V of a double helix of N base pairs is thus approximated by the multi-dimensional parabola,

(3)

where _i is the instantaneous fluctuation of the ith "step" parameter from its equilibrium value and f_ij is an element of a 6(N  1) × 6(N  1) matrix of elastic force constants, F, i.e., the terms describing the potential deformability of the six parameters in each of the N  1 base-pair steps.

The kinetic energy K is similarly expressed in quadratic form in terms of _i, the first derivative of _i with respect to time, and the weighted "mass" coefficients h_ij of the 6(N  1) × 6(N  1) matrix H (defined below):

(4)

If the _i and _i are collected respectively in the vectors  and , the equations of motion simplify to:

(5)

with a periodic solution of the form:

(6)

Here _n is the frequency, _n the phase angle, _n the amplitude of the nth normal mode, and A_in the fractional contribution to the ith "step" parameter from the nth normal mode. (In practice the components of and the frequencies _n are obtained by introducing the vector of normal coordinates Q defined by  = AQ, where the elements of Q are _n cos(_nt + _n) and A is a 6(N  1) × 6(N  1) transformation matrix. Eq. 5 can then be rewritten as the eigenvalue expression, HA = FA, with A^THA equated to the identity matrix of order 6(N  1) and the normal mode frequencies and directions found from the diagonalization of F, i.e.,  = A^TFA, where A^T denotes the transpose of A and  is a diagonal matrix with elements _nn = . The phase angles, however, cannot be determined with this procedure.) The fluctuations of each "step" parameter _i at time t are thereby expressed as a linear combination of harmonic oscillators, the energies of which are proportional to . The contribution of each mode to the variation of individual "step" parameters decreases rapidly with increase in _n, so that relatively few (low frequency) modes are responsible for the large-scale deformations of the molecule.

Kinetic energy

The kinetic energy coefficients in Eq. 4 incorporate the mass m_a and the Cartesian coordinates r_a of individual atoms in the DNA through the relationship:

(7)

To evaluate the h_ij, we take advantage of analytical expressions for the (r_a/_i) developed to treat the normal mode dynamics of a system of two molecules, each of which moves in dihedral angle space (Braun et al., 1984; Higo et al., 1985). The set of rigid body parameters used here to relate neighboring base-pair planes is identical in form to the variables previously used to describe the relative global positions and orientations of different molecules. The internal atomic motions of DNA are effectively separated from the overall rotations and translations of the molecule by expressing the atomic displacements in terms of an embedded coordinate frame, chosen according to the so-called Eckart condition (Eckart, 1935) to minimize the mass-weighted square atomic displacement of DNA (McLachlan, 1979) in its equilibrium rest state and in a state where one of the "step" parameters is altered from its minimum energy value (Noguti and G, 1983).

Because the atomic motions r_a brought about by "step" parameter change _i are assumed to be small in this common frame (Eckart, 1935), the kinetic energy term (Eq. 4) effectively restricts the length of DNA which can be studied by normal mode analysis. This limitation does not apply to the mesoscopic chains lengths considered below, where typical room temperature fluctuations of individual base-pair steps, such as the ±5° perturbations of Roll which raise the potential energy by ~k_BT/2, limit atomic movement to 0.02-0.07 persistence lengths. The persistence length of mixed sequence DNA under ambient aqueous salt conditions, by comparison, is the same magnitude as the contour length of a 150 to 200-bp chain (Hagerman, 1988; Smith et al., 1992, 1996; Bustamante et al., 1994; Bednar et al., 1995; Baumann et al., 1997). The limitation on chain length is further discussed below in the analysis of DNA end-to-end dimensions.

Both backbone and base atoms must be included in the evaluation of the h_ij. The sugar-phosphate backbone is incorporated in the present calculations by the superposition of a canonical B-form 5'-nucleotide helical fragment (Chandrasekaran and Arnott, 1989) in the reference frame of each base. Because complementary base-pair parameters are fixed at ideal planar values, each complementary nucleotide pair is thereby treated as a rigid body and the small variations in backbone conformation that accompany fluctuations in the geometry of neighboring base pairs are ignored.

Force field

The potential energy coefficients in Eq. 3 are taken from knowledge-based elastic functions extracted from the three-dimensional arrangements of DNA base-pair steps in protein-DNA crystal structures (Olson et al., 1998):

(8)

The g_ij in this expression are the force constants impeding deformations of the XpY dimer, and the _i correspond to the instantaneous fluctuations of each of the six "step" parameters of that dimer from the B-DNA rest state, i.e., ₁ = Tilt_XY, ₂ = Roll_XY, ₃ = Twist_XY, ₄ = Shift_XY, ₅ = Slide_XY, ₆ = Rise_XY (full details of the potential function are available in Table S-1 in the Supplementary Material). If these elements are collected respectively in the 6 × 6 force constant matrix G_XY and the 6 × 1 vector _XY, the dimer deformation energy V_XY reduces to (1/2)_XY^TG_XYXY. The choice of translational parameters and the aforementioned approximation of T in Eq. 2 by the constant matrix T_o introduce a linear relationship between the _XY and the corresponding elements of the longer 6(N  1) vector used in Eq. 5. That is, _XY = C_XYXY, where the molecular deformation vector (Eq. 5) is built up from the (N  1) sets of constituent _XY dimer fluctuations and C_XY is a 6 × 6 constant matrix characteristic of the mean orientation of the XY dimer in B-form DNA. The potential energy of the DNA as a whole can be expressed in the form of Eq. 3 by constructing the 6(N  1) pseudodiagonal matrix of force constants F from the F_XY = C_XY^TG_XYC_XY corresponding to the DNA base-pair sequence. For example, the f_ij of the self-complementary (CGTACG)₂ hexamer duplex would be collected in the 30 × 30 array,

(9)

where the boldface 0 values in this expression are 6 × 6 null matrices.

Because the lowest energy conformation of a given DNA sequence is self-evident from the knowledge-based force field, i.e., the minimum energy three-dimensional structure is defined by the set of average "step" parameters, the time-consuming energy minimization step normally carried out before normal mode analysis does not need to be performed in the present calculations.

The potential functions introduced below are first approximations of the sequence-dependent structure of DNA based on the observed conformations of selected dimers in high resolution protein-bound and B-DNA crystal structures. Although mean (equilibrium) sequence-dependent parameters have not significantly changed since first reports (Gorin et al., 1995; Olson et al., 1998), the local force constants are expected to change as new structural data accumulate. Moreover, there are not yet sufficient crystallographic data to address DNA deformability beyond the dimer level despite suggestions (Nadeau and Crothers, 1989) of longer-range organization of DNA structure.

End-to-end dimensions

The unperturbed, mean-square end-to-end distance of DNA, r²₀, is obtained by three independent approaches. First, exact values of r²₀ are computed using standard matrix methods for the mean extension of a DNA of N base pairs, with r²₀ equated to the average scalar product of the end-to-end vector, r · r₀, and r defined as

(10)

Here T_k,k+1 is the transformation matrix that relates coordinate frames in successive base-pair planes (Eq. 1) and v_k is the virtual bond vector that connects base-pair centers (Eq. 2). Boltzmann averages of the T and v at each base-pair step are incorporated in dimeric generator matrices, the serial product of which yields all possible configurations of the selected chain sequence. The reader is referred elsewhere for the precise formulation of these expressions (Flory, 1969; Maroun and Olson, 1988; Marky and Olson, 1994). The zero subscript on r²₀ and other (unperturbed) averages denotes the omission of chain self-intersection and solvent-polymer interactions in their calculation (Flory, 1969). The DNAs studied below are sufficiently short and stiff so that excluded volume effects are automatically avoided.

Second, we express r²₀ in terms of the fluctuations r of the end-to-end vector with respect to the position r^o of terminal base pairs in the (minimum energy) rest state,

(11)

Here the superscript ^o denotes (constant) rest state values. The value of r, which is determined by the fluctuations of "step" parameters, is written as a Taylor series expansion in the _i,

(12)

If the series is approximated by the first term, the average displacement vector r₀ is zero because _i = 0 and the mean-square end-to-end distance is greater than the separation of chain ends in the rest state (because the self product, r · r₀, in Eq. 11 must be positive). Because such a result is inconsistent with the end-to-end displacement of a naturally straight, inextensible DNA, we consider up to second-order terms of the Taylor series (Eq. 12) in the calculation with the result,

(13)

The first and second derivatives in these expressions are evaluated numerically. The values of _ij are taken directly from the experimental data used to derive the knowledge-based force field (Olson et al., 1998). Such averages would have to be determined numerically over all normal modes if an all-atom force field were used. It should also be noted that the first derivatives that appear in Eq. 13, (r/_i), differ from similar terms employed in Eq. 7 to compute the kinetic energy. As explained above, the latter derivatives are computed in a carefully chosen coordinate frame that essentially eliminates rotational and translation motions of the molecule. The variation of chain ends required to evaluate Eq. 13, by contrast, is performed in a common global coordinate frame embedded in the plane of the first base pair.

Finally, we compute the root-mean-square end-to-end distance of DNA molecules which are confined to the set of spatial configurations associated with one or more normal modes. That is, we use the set of "step" parameters associated with mode n at time t (Eq. 6) to generate the end-to-end vector r(n, t) and the Boltzmann factor (n, t) of each of the polymeric states that comprise certain global motions and then evaluate r²₀ as an energy weighted sum over representative configurations:

(14)

In practice, we set an upper limit of 7.5 k_BT on the potential energy, V(_i(n, t)) (Eq. 3), thereby excluding configurations where (n, t) = exp[V/k_BT]  exp[7.5].

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS GLOSSARY REFERENCES

Elastic features of DNA repeating polymers

We start by determining the normal modes of four simple repeating polymersthe poly dA · poly dT and poly dG · poly dC homopolymers with a monomeric repeating unit (the AA · TT or GG · CC dimer) and the poly d(AT) and poly d(GC) alternating copolymers made up of sequential purine-pyrimidine (AT · AT or GC · GC) and pyrimidine-purine (TA · TA or CG · CG) base-pair steps. We classify the motions, after Matsumoto and G (1999), on the basis of the global structural changes revealed through computer visualization of the eigenvectors associated with each mode (i.e., Eq. 6). The color-coded spectrum of lowest frequency bending, twisting, and stretching modes of a 200 bp chain is presented in Fig. 1 for each of the four polymers and schematic illustrations of some of these motions are given in Fig. 2. Numerical values of the frequencies are tabulated in the Supplementary Material (Table S-2). We focus on the lowest frequency modes because the amplitudes of the normal modes and their consequent importance to the overall motions of DNA diminish at higher frequency.

View larger version (26K): [in this window] [in a new window]   FIGURE 1   Color-coded spectrum of lowest frequency bending (blue), twisting (green), and stretching (red) modes of 200-bp segments of the poly dA · poly dT and poly dG · poly dC homopolymers and the poly d(AT) and poly d(GC) alternating copolymers. As the frequency increases, the amplitude of each normal mode and its consequent importance to the overall motions of the DNA becomes lower.

View larger version (18K): [in this window] [in a new window]   FIGURE 2   Schematic illustration of representative low frequency normal modes of an elastic rod. The arrows point in the directions of bending, twisting, and stretching motions in each mode.

It should be noted that these low frequency modes would not be observed in a real system due to viscous effects, which we do not consider here. While our method does not necessarily describe the time-development of DNA chains in a real system, we show below that these low normal mode frequencies are closely related to the global dynamic properties of DNA, such as the bending rigidity, twisting rigidity, and stretching rigidity. Thus, the low normal mode frequencies presented here are useful indicators of the dynamic properties of DNA.

Each of the pure (planar) bending modes (blue lines at the left of Fig. 1) is doubly degenerate, describing mutually perpendicular motions with equivalent (bending) energy. The superposition of these modes leads to the global isotropic bending characteristic of an ideal elastic rod. In contrast to the bending motions, the twisting modes (green lines) are higher in frequency, more widely spaced, and more sensitive to polymer composition. Specifically, the frequencies of the twisting modes of the simulated alternating copolymers and the poly dG · poly dC homopolymer are lower in value than those of the poly dA · poly dT model. The stretching modes (red lines), which are also sensitive to base-pair sequence, tend to be even more widely spaced and to occur at still higher frequencies than the twisting modes.

The rod-like behavior of poly dA · poly dT is further illustrated in Fig. 3 where the chain length dependence of the computed low frequency modes is superimposed on the predicted variation (Strutt, 1894) of the bending, twisting, and stretching frequencies of a homogeneous elastic rod. (The frequency of the nth bending mode of an elastic rod of length L with uniform circular cross-section is given by  = C_bp, where C_b is a material constant and p_n satisfies the condition p_nL = (4.730, 7.853, 10.9965, ···) for (n = 1, 2, 3, ...) (Strutt, 1894). The twisting and stretching frequencies follow a similar form, varying directly with n and inversely with L, i.e.,  = (n/2L)C_t and  = (n/2L)C_s.) The frequencies of the poly dA · poly dT motions (scatter points) are normalized with respect to the lowest frequency bending, twisting, and stretching modes of a 200-bp chain. The normal mode frequencies of the elastic rod (smooth curves) are similarly expressed as multiples of the lowest frequency motions of an ideal rod with the same contour length (199 base-pair steps × 3.27 Å per step = 650 Å). The agreement between computed and theoretical values is remarkably close, with only minor discrepancies in the relative bending frequencies at short chain lengths where the theory is known to break down. Similar agreement between the base-pair level representation of DNA motions and elastic rod behavior is obtained for the low frequencies modes of poly dG · poly dC and the two alternating copolymers (data not shown).

View larger version (22K): [in this window] [in a new window]   FIGURE 3   Superposition of the computed dependence of representative low frequency bending, twisting, stretching modes of poly dA · poly dT (scatter points) on chain length and the predicted variation (thin solid lines for the lowest frequencies, dashed lines for the second lowest frequencies, and thick solid lines for the third lowest frequencies) for an elastic rod with the same contour length. Frequencies are normalized with respect to the lowest frequency of each type of motion in a 200 bp-chain.

The mechanical constants associated with the overall elastic behavior of the four repeating polymers are collected in Table 1. These values are obtained by substituting the computed normal mode frequencies in the classical expressions for the normal modes of an ideal elastic rod (Strutt, 1894). Specifically, the bending constant A, the twisting constant C, and the stretching constant Y are given by

(15)

where the _n are the computed frequencies of the nth bending (b), twisting (t), or stretching (s) mode, M is the total mass of the molecule, L is the length of the DNA helical axis (see Kosikov et al., 1999, and Matsumoto and G, 1999 for computational details) and I_M is the moment of inertia per unit length around the twisting axis. As evident from the table, the computed values of the bending constant A and the persistence length a derived from A on the basis of elastic rod theory (a = A/k_BT with T = 298 K) are approximately independent of sequence and in excellent agreement with values extracted from experimental studies of mixed sequence DNA (Hagerman, 1988; Smith et al., 1992, 1996; Bustamante et al., 1994; Bednar et al., 1995; Baumann et al., 1997). The elastic constants, however, may change as the local force constants are refined.

The sensitivity of the calculated twisting constants to base-pair sequence may account in part for the wide range of experimentally derived values of C (Barkley and Zimm, 1979; Thomas et al., 1980; Hurley et al., 1982; Millar et al., 1982; Horowitz and Wang, 1984). The differences in measured twisting constants, usually attributed to experimental effects, such as the bending strain that accompanies loop closure (Heath et al., 1996), may also reflect the intrinsic deformability of the constituent base pairs. The greater overall twisting stiffness of poly dA · poly dT compared to the other polymers in Table 1 stems in part from the substantially higher cost of distorting AA dimers compared to other base-pair steps in the knowledge-based potential (see Table S-1 in the Supplementary Material and discussion below).

The stretching constants Y extracted here, however, are two- to threefold greater than values deduced from recent single-molecule manipulations of DNA (Smith et al., 1996; Baumann et al., 1997; Wang et al., 1997; Bouchiat and Mezard, 1998). As discussed below, the relative ease of global stretching reflects a number of local factors, in addition to the force constants that govern the axial, i.e., van der Waals', separation of neighboring base-pair planes. In other words, the global stretching of DNA is not simply a function of Rise.

Matsumoto and G (1999) previously derived DNA elastic constants on the basis of a normal mode analysis in the dihedral angle space of short (24-36 bp) double helical molecules, with poorer correspondence between computation and experiment than the present work, i.e., the previously computed bending and twisting constants are respectively three to ten times greater than the corresponding values in Table 1. The improved agreement found here seemingly reflects our use of knowledge-based potential energies rather than an atomic force field. Optimization of DNA conformation on the basis of the detailed interactions of all atoms typically yields a large number of closely spaced minimum energy substates of similar structural character (Poncin et al., 1992). The normal mode analysis of DNA in dihedral angle (and Cartesian) space is performed with respect to one of these minimum energy states and the transitions between closely related, competing minima are not considered. The statistical approach used to generate our elastic force field, by contrast, reflects a large number of known crystal structures, each of which is analogous to one of the minimum energy substates found through all-atom energy minimization. Therefore, any fluctuation on our elastic energy surface is comparable to a transition between different minimum energy points derived with an atomic force field. The correspondence between the mechanical constants in Table 1 and values extracted from experiment suggests that transitions between conformational substates have a significant influence on the global properties of DNA.

A similar argument is made by Song and Schurr (1990) in accounting for the unusually large values of DNA persistence length (~2100 Å) deduced from measurements of transient electric dichroism. They offer three possible reasons for the large discrepancy between their results and the persistence lengths for mixed sequence DNA (~500 Å) obtained with other experimental approaches. One explanation, which they favor, is that the potential energy of DNA bending is not a smooth quadratic function and that the energy exhibits several discrete minima separated by barriers. In the short-time scale of their experiment the DNA is trapped in one of the minima and the persistence length thereby becomes larger.

It should be noted that there is an open controversy over the ionic strength dependence of the persistence length and the lower limit of the persistence length of mixed sequence DNA at ambient salt conditions. We cite a DNA persistence length of 500 Å in Table 1 consistent with a variety of independent measurements. As noted above, larger values of the persistence length are deduced in some work, e.g., (Song and Schurr, 1990). There is similar disagreement among elastic constants obtained for specific DNA sequences (see Table S-3 in the Supplementary Material) which may be related in part to the discussion above, i.e., chains may sample only a small part of the energy surface in some experiments but may explore vast areas of conformation space in others.

Sequence effects on global motion

Figs. 4-6 illustrate the deformations of local "step" parameters, which are collectively responsible for selected low frequency normal modes of the aforementioned repeating polymers. These plots, for DNA chains of 120 bp, capture the sequential fluctuations of each conformational variable at the moment when the potential energy of the molecule is raised by k_BT/2; the fluctuations of all parameters are reversed a half cycle later of the mode. As detailed below, the mechanisms used to effect overall chain motions reflect the conformational properties of the constituent dimer units.

View larger version (42K): [in this window] [in a new window]   FIGURE 4   Deformations of local angular "step" parameters, which are collectively responsible for the lowest frequency in-plane bending motions of a 120-bp poly dA · poly dT duplex: (top, middle) degenerate modes of lowest frequency, nb = nb' = 1; (bottom) one of the two second lowest frequency modes, nb = 2. Plots illustrate the sequential fluctuations of Tilt (thin solid line), Roll (dashed line), and Twist (thick solid line) at the moment when the potential energy of the molecule is raised by kBT/2; fluctuations are reversed a half cycle later of the mode.

View larger version (33K): [in this window] [in a new window]   FIGURE 5   Fluctuations of local angular "step" parameters, which are collectively responsible for the lowest frequency global twisting modes (nt = 1) of 120-bp fragments of poly dA · poly dT, poly d(AT), poly dG · poly dC, and poly d(GC). See legend to Fig. 4.

View larger version (36K): [in this window] [in a new window]   FIGURE 6   Fluctuations of local translational "step" parameters associated with the lowest frequency global stretching (ns = 1) of 120-bp fragments of poly dA · poly dT, poly d(AT), poly dG · poly dC, and poly d(GC). Plots illustrate the sequential fluctuations of Shift (thin solid line), Slide (dashed line), and Rise (thick solid line) at the moment when the potential energy of the molecule is raised by kBT/2. See legend to Fig. 4.

Bending

1

1

Twisting

Stretching