Ultraviolet (UV) difference spectroscopy

Figure 1 shows the UV difference spectra of T7 DNA at each of the several temperatures relative to a matching sample at 3 °C (Sarocchi and Guschlbauer, Reference Sarocchi and Guschlbauer1973) (SG). These results were ascribed to the broad pre-melting transition of DNA. It was suggested originally that the two states responsible for the isosbestic points in the UV difference spectra corresponded to the unprotonated and protonated forms of oligodC⋅oligodG subsequences, but that suggestion was abandoned after further studies showed that polydC⋅polydG remained unprotonated at neutral or higher pH, where the measurements were made. The existence of isosbestic points at Δɛ = 0 in these difference spectra near 240 and 265 nm are strong evidence for a shift of population between two different states that exhibit a T-independent difference between their intrinsic molar absorbances at any given wavelength, as T is gradually increased from 3 to 45 °C, and possibly even higher to 60 °C at the longer wavelengths. The non-conforming behavior of the 60 °C sample at wavelengths below 250 nm might reflect some variation of the difference in the intrinsic short-wavelength absorbance between the two states with increasing temperature above 45 °C. The modest displacement of sample 2 from the isosbestic point near 263 nm is likely due to incomplete equilibration upon increasing the temperature from 3 °C. Abundant evidence that the transition is very slow in the low-T range between 3 and 22 °C was presented. The observed absorbance differences are quite small, amounting to ~1% of the maximum absorbance near 260 nm. Such changes in UV absorbance and also CD spectra were entirely reversible.

Fig. 1. UV difference spectra for T7 DNA in 0.1 M NaCl, pH 7.6, at various temperatures. Δɛ is the difference in molar absorbance between the sample at temperature T and a matching sample held at 3 °C. Data were recorded using a dual beam instrument with a reference sample held at 3 °C, while T of the other sample was raised in a uniformly increasing manner. Curves 1–6 correspond to temperatures: 3, 22, 35, 45, and 60 °C, respectively. This is Fig. 1B of Sarocchi and Guschlbauer (Reference Sarocchi and Guschlbauer1973).

X-ray diffraction from crystals

A high-resolution (0.96 Å) electron density map for crystals of the duplex decamer, CCAGGCCTGG, which previously had been examined at 1.6 Å resolution and refined to a single structure, was analyzed (Maehigashi et al., Reference Maehigashi, Hsiao, Woods, Moulaei, Hud and Williams2012). The new data and methods of the analysis revealed extensive heterogeneity of DNA conformation, in particular, the existence of multiple DNA structures coexisting at equivalent positions in different unit cells of the same crystal with two distinct structures being more or less dominant. Structural differences between the two dominant states were greater for the backbones than for most of the bases. Heterogeneity in the positions and occupancies of water molecules and ions was also reported. As stated by the authors, these ‘high-resolution data reveal that intermediate and low-resolution data can give a false sense of homogeneity.’ In light of that, structural heterogeneity of a given DNA in the same crystal might be far more widespread than is commonly supposed. It is unlikely to be less in solution.

X-ray scattering from DNAs labeled with gold colloids in solution

A group at Stanford presented the results of three different small-angle X-ray scattering (SAXS) studies of solutions containing duplex DNAs of various lengths from N = 10 to 35 bp that were site-specifically labeled at their ends and other positions with gold colloids (Mathew-Fenn et al., Reference Mathew-Fenn, Das and Harbury2008, correction, 2009; Shi et al., Reference Shi, Herschlag and Harbury2013; Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016). All three studies reported that the variance of the total end-colloid to end-colloid distance contained a small constant part plus a part proportional to N ². The most recent and precise study (anomalous SAXS) yielded an average rise per bp, h ₀ = 0.323 ± 0.01 nm, and concluded that each bp contributed about 0.025 nm (in-phase) to the total standard deviation of the end-to-end distance distribution (Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016). This is consistent with a largely all-or-none variation in length, as would be expected for a highly cooperative two-state equilibrium, wherein the two states are roughly equally populated and their rises per bp are: (0.323 + 0.025) = 0.348 nm and (0.323 – 0.025) = 0.298 nm.

The conclusions of the original paper (Mathew-Fenn et al., Reference Mathew-Fenn, Das and Harbury2008) were quickly challenged by theorists (Mazur, Reference Mazur2009; Zheng et al., Reference Zheng, Czapla, Srinivasan and Olson2010), who found a component of the colloid-to-colloid variance that varied as N ² for models that underwent bending but did not exhibit a cooperative stretching transition. This quadratic behavior arose from end-to-end fluctuations caused by bending, which were significantly amplified by the vectors between the ends of the DNA and the mean positions of their colloid centers, and further amplified by variations in colloid position about its mean position. A convincing rebuttal that included many additional scattering data for gold colloids attached at different positions on the DNA, and analyses using models that incorporated twisting as well as bending were presented (Shi et al., Reference Shi, Herschlag and Harbury2013). The constant contribution to the variance of the colloid-to-colloid distances was very small, indicating that the distribution of positions of the gold colloids with respect to their points of attachment was far narrower than those assumed by Mazur (Reference Mazur2009) and Zheng et al. (Reference Zheng, Czapla, Srinivasan and Olson2010), thereby ruling out those explanations. The ‘knowledge-based’ model used by Zheng et al. (Reference Zheng, Czapla, Srinivasan and Olson2010) gave a poor fit to both the original and newer data. Predicted periodic oscillations of the variance with N due to the displacement of the colloids from the helix-axis (Becker and Everaers, Reference Becker and Everaers2009) were observed. Even after including effects of bending, twisting, and offset of the colloids from the helix-axis, there remained an excess variance that increased as N ², and a two-state all-or-none model still gave the best fit. A subsequent study using improved methods and smaller colloids yielded a significantly larger N ² contribution to the variance than that of the first study (Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016).

A weakly cooperative two-state model proposed earlier (Storm and Nelson, Reference Storm and Nelson2003) was adapted to treat numerically the properties of short DNA chains (Padinhateeri and Menon, Reference Padinhateeri and Menon2013), and was claimed to account for the X-ray scattering data of Mathew-Fenn et al. (Reference Mathew-Fenn, Das and Harbury2008). However, it failed to do so in multiple ways. The distribution functions were calculated for the end-to-end distances of the DNA alone and compared with the distribution of measured colloid-to-colloid distances. The end-to-end distances do not include the vector from each end of the DNA to the center of its gold colloid and the variation allowed by its tether, which contributes to the peak position and width of the X-ray scattering profile. Moreover, the model cannot fit the experimental distance distributions for different lengths without also decreasing the statistical weight of the less populous state by more than 2-fold as N increases from 30 to 35 bp. The calculated variance of the end-to-end distance varies as ~N ^3.5 instead of the ~N ² observed.

Nuclear magnetic resonance (NMR) of small DNAs in solution

NMR evidence for a dynamic two-state equilibrium between two different B-form structures in a 4 bp DNA was reported (Reid et al., Reference Reid, Salisbury, Brown and Williams1985). Later evidence for a two-state dynamic equilibrium in a 7 bp heteroduplex DNA was reported, although no detailed structural information was obtained (Manoharan et al., Reference Manoharan, Gerlt, Wilde, Withka and Bolton1987). Slow conformational dynamics of the adenines in T–A steps with all possible nearest neighbors in duplex dodecamers was also reported, but detailed structural information was lacking (McAteer et al., Reference McAteer, Ellis and Kennedy2000). Slow site-exchange dynamics (t ≥ 30 ms) of the thymine complementary to the adenine at a TpA step of a 16 bp duplex DNA provided further evidence of a slow structural change (Diaz, Reference Diaz2002). The paucity of NMR solution structures pertaining to different coexisting structures in the same duplex DNA is likely due at least in part to the selection of DNAs with sufficiently narrow and well-separated cross-peaks in the two-dimensional (2D) nuclear overhauser effect spectroscopy and correlated spectroscopy spectra to enable structure determination. In the lab of the late Prof. Brian Reid, for every short (6–16 bp) duplex oligomer, whose full 3D structure was determined, another ~20+ such oligomers were synthesized, but their 2D NMR spectra were judged unsuitable for structure determination. In any case, the majority of sequences studied anywhere by NMR were capped at their ends with multiple G⋅C bp to provide stability against dissociation and fraying, and that may have tilted the prevailing equilibrium heavily in favor of one or another particular conformation.

Formulation of the model

The statistical weight of a bp in the a state, when preceded by a bp in the a state or by solvent, is 1.0. However, the statistical weight of a bp in the a state, when preceded by a bp in the b state, is J × (1.0) = J, where J is the cooperativity parameter. J is expected to be <1.0, reflecting an extra free energy cost associated with the formation of a junction between the b and a conformations. Similarly, when a bp in the b state is preceded by a bp in the b state or by solvent, its statistical weight is B, and when preceded by a bp in the a state, its statistical weight is J × B = JB. For example, the statistical weight of DNA containing 10 bp in the sequence of states, aaabbbbaaa, is 1³JB ⁴J1³. The quantity B is the equilibrium constant for an a⇔b transition of a single bp (without forming a new junction). The statistical weight of any given sequence of states for an N bp DNA is a product of N factors, each of which is either 1, J, B, or JB. The total unnormalized statistical weight of the model is the sum over the statistical weights of all possible sequences and is called the partition function, χ. Using the rules of matrix multiplication and theory of linear transformations, χ can be expressed as

(1)$$\chi \, = \,\lsqb {1,B} \rsqb ({\bi M}^{N-1} )\,\left[ {\matrix{ 1 \cr 1 \cr } } \right] = \,\lsqb {1,B} \rsqb {\bi S}\left\vert {\matrix{ {\lambda_ +^{N-1} } & 0 \cr 0 & {\lambda_-^{N-1} } \cr } } \right\vert {\bi S}^{-1}\left[ {\matrix{ 1 \cr 1 \cr } } \right]$$

where the so-called transfer matrix is given by

(2)$${\bi M}\equiv \left\vert {\matrix{ 1 & {JB} \cr J & B \cr } } \right\vert $$

the λ _± = (1/2) (B + 1 ± ((B − 1)² + 4BJ ²)^1/2) are eigenvalues of M, and S is a similarity transformation matrix that diagonalizes M. That is, S⁻¹ MS = λ, wherein λ is the diagonal matrix with elements λ ₊ and λ ₋. The columns of S are right eigenvectors of M. The normalized statistical weight of any configuration of the a and b states in the entire DNA is just the unnormalized weight of that configuration divided by χ. The average number of bp in the b state is readily shown to be, 〈n_b〉 = (B/χ) × (dχ/dB), and the average number of junctions is, 〈n_J〉 = (J/χ) × (dχ/dJ), both of which are readily calculated. When N is sufficiently large that λ₋^{N −1} is negligibly small compared to λ ₊^{N −1}, the average fraction of bp in the b state is readily found to be

(3)$$ {\,f_b} \equiv \left\langle {{n_b}} \right\rangle /N\, = \,\left( {1/2} \right)(1\, + \,\left( {B\, - \,1} \right)/{({\left( {B\, - \,1} \right)^2}\, + \,4B{J^2})^{1/2}})$$

For arbitrarily small values of N, such that λ₋^{N −1} is no longer negligible, the quantities, χ, 〈n_b〉 and f_b can be expressed algebraically and evaluated numerically using Mathematica. Average properties, such as length and inverse torsion, and bending elastic constants are readily calculated from f_b and the intrinsic values of such properties that are associated with the individual a and b states. The individual a and b states are assumed to be non-stretchable but have different rises per bp, so any apparent stretching of the DNA results from shifting the equilibrium constant and the relative populations of the two states.

If the DNA is held under tension, F, then the average additional free energy change between the longer b state and the shorter a state of any given bp is approximately, w = − Fδ〈cos θ〉, where δ is the difference in length between the b and a states, and θ is the angle between the (helix) axis of the particular bp and the applied force. The average of cos θ is taken over all bps in the DNA. The equilibrium constant in the presence of the force is approximately

(4)$$B_F = B_0\exp \lsqb + F\delta \left\langle {\cos \theta } \right\rangle /k_{\rm B}T\rsqb $$

where B ₀ is the equilibrium constant in the absence of force and k _BT is the thermal energy. For the values of F and δ considered here, the positive quantity in the exponent is far <1.0, but nevertheless has a substantial effect due to the large cooperativity (very small value of J), as will be seen.

Estimation of model parameters from experiment

Several remarkable theoretical and experimental advances enabled quantitative estimation of the model parameters. First came an extremely accurate theory of relative extension versus force for a non-stretchable worm-like coil (WLC) with a fixed contour length and bending rigidity (Marko and Siggia, Reference Marko and Siggia1995). The primary result was a numerical algorithm for extremely precise calculations. Then came analytical expressions for both relative extension versus force and torque versus force for a slightly twisted, non-stretchable but twistable, WLC (TWLC) model (Moroz and Nelson, Reference Moroz and Nelson1997, Reference Moroz and Nelson1998). These expressions were accurate for small torques and sufficiently high forces. For untwisted DNAs at forces, F ≥ 0.25 pN, the Moroz–Nelson (MN) formula for relative extension versus force is very accurate (<0.2% deviation from Marko–Siggia theory; Schurr, Reference Schurr2015). The MN analytical torque expression agrees within statistical errors (4%) with simulations of a twisted discrete model of a non-stretchable WLC model at various forces from 0.25 to 4.0 pN (Chou et al., Reference Chou, Lipfert and Das2014).

Throughout this paper, the term ‘stretching’ means lengthening of the helix-axis, or increasing the rise per bp, whereas the term ‘extension’ means increasing the average projection of the helix-axis onto the direction of the pulling force.

Although the MN formula for relative extension of an untwisted DNA versus force was derived for a model with constant contour length L and constant persistence length P, it should also apply at least approximately to a stretchable model, wherein both the contour length and the persistence length vary with force. It takes the form

(5)$$ z_F/L_F\, = \,(1\,-\,\lpar {0.5} \rpar /(P_FF/k_{\rm B}T\,-\,1/32)^{1/2})$$

where z_F, L_F, and P_F ≡ h_F κ_βF/k _BT, are, respectively, the extension, contour length, and persistence length, and h_F and κ_βF are, respectively, the mean rise per bp and bending elastic constant of the springs between bps, of a DNA held at force F.

A decade later came two extraordinary experimental advances. First was the acquisition of precise relative extension versus twist data for DNAs held at many different forces and twists in the low-force regime from 0.18 to 3.9 pN, from which the effective torsional rigidities were extracted by a thermodynamic path integration method (Mosconi et al., Reference Mosconi, Allemand, Bensimon and Croquette2009). The second extraordinary advance was a direct measurement of the torque on a single DNA molecule held at various tensions from 0.2 to 6.5 pN as a function of twist in the small-twist limit by two different magnetic tweezer-based methods (Lipfert et al., Reference Lipfert, Kerssemakers, Jager and Dekker2010, Reference Lipfert, Wiggin, Kersemakers, Pedaci and Dekker2011). The effective torsional rigidity was obtained directly from the slope of the torque with respect to the overall twist angle in the small-twist limit.

Unfortunately, the contour length of the DNA in most, if not all, single-molecule measurements on DNAs is not precisely known. The effective contour length, z _max, is typically estimated by fitting an analytical formula (Bouchiat et al., Reference Bouchiat, Wang, Allemand, Strick, Block and Croquette1999), which closely approximates the Marko–Siggia theory to the experimental z_F versus F data to determine best-fit values of P and z _max, from which z_F/z _max is reckoned.

The parameters of the two-state model were estimated by fitting the model to various data while invoking assumptions that were later found to be modestly incorrect. For example, in determining the variation of L_F with increasing force, the persistence length of DNA in 0.1 M NaCl was originally assumed to remain a constant, P = 50 nm, independent of the force. Later work indicated that P_F actually varies significantly, though modestly, with increasing force. The various fitting protocols were described and implemented in previous studies (Schurr, Reference Schurr2015, Reference Schurr2017, Reference Schurr2019a), and the resulting parameters obtained are listed in Table 1.

Table 1. Fixed parameters of the cooperative two-state model

The value of δ is close to 0.05 nm estimated from X-ray scattering studies of short DNAs labeled with gold colloids (Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016). The value of h ₀ yields, h_av = 0.323 nm, for an unstrained DNA in 0.1 M NaCl at 293 K, identical to the average value, 0.323 ± 0.01 nm, determined for short DNAs under such conditions (Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016).

The α_a and κ_βa are effective Hooke's Law torsion and bending elastic constants, respectively, of the effective springs between bps in the a conformation, and α_b and κ_βb are the corresponding constants for springs between bps in the b conformation. These elastic constants apply to the twisting and bending angles that orient a coordinate frame of one subunit in the coordinate frame of its preceding subunit. In lieu of a full discussion in terms of the frame-to-frame Euler rotations (Wilcoxon and Schurr, Reference Wilcoxon and Schurr1983, Appendix), the following simplified description is offered. Imagine a straight filament comprising identical stacked cylindrical disks, each containing the atoms of a single bp step and centered on the common straight symmetry axis. If a given disk is rotated around its symmetry axis by an angle ϕ with respect to its preceding neighbor, then ϕ (radians) is the twisting angle of the torsion spring between those two disks. The twist required to attain thermal energy in a particular torsion spring is ϕ = (2k _BT/α)^1/2, which at T = 293 K yields ϕ  = 0.135 rad (7.7°) for α_a and 0.0810 rad (4.6°) for α_b.

For a DNA that is not straight, the bending angle β between two adjacent disks is the angle between their symmetry axes. The apparent bending elastic constant, $\kappa _\beta$, of DNA consists of fast-relaxing, slowly-relaxing, and non-relaxing parts (Naimushin et al., Reference Naimushin, Fujimoto and Schurr2000; Schurr, Reference Schurr2019a). For DNA in 0.1 M NaCl at 293 K, the fast-relaxing part of the persistence length is, P_f = 170 nm (Okonogi et al., Reference Okonogi, Reese, Alley, Hopkins and Robinson1999), which is assumed to apply to both a and b states. The corresponding fast-relaxing part of the bending elastic constant of the a state is given by κ _βaf = (_f/h_a)k _BT and a corresponding expression applies for the b state. For these elastic constants, the bending angle β required to attain thermal energy at 293 K is β _a = (2k _BT/κ _βaf)^1/2 = 0.059 rad (3.4°) for the a state and β _b = 0.064 rad (3.7°) for the b state.

The fast-relaxing part of the bending rigidity of DNA, A_f = (k _BT)P_f = 6.88 × 10^–28 Nm², is the appropriate quantity for comparison with straight cylindrical rods of other materials. The radius of a cylindrical sheath that just encloses the outermost atoms of a straight DNA is R = 1.0 nm. Even though a significant fraction of the volume inside that cylinder is water and ions, the entire cylinder will be regarded as ‘DNA’ for the purpose of this discussion. The bending rigidities of cylinders of the same radius, R = 1.0 nm, are obtained for other materials according to A = E (π/4) R ⁴, where E is the known flexural (Young) modulus of the material. The bending rigidity of the DNA ‘cylinder’ of 1.0 nm radius is found to be 1.6-fold greater than that of high-density polyethylene, 0.8-fold smaller than that of polypropylene, 0.34-fold smaller than that cellulose acetate, and 0.29-fold smaller than that of solid nylon. The metals, magnesium, aluminum, and iron are, respectively, 51-, 79-, and 228-fold stiffer than the DNA cylinder. If the other materials had the same shape as DNA in its cylinder, their bending rigidities would be decreased substantially, perhaps by 2-fold or more, in which case DNA would appear more like cellulose acetate or solid nylon. In regard to rapidly relaxed bending, DNA lies in the range of ordinary polymeric materials.

The mean domain size at the midpoint of the transition (B = 1.0) is 1 + 1/J = 224 bp. The equilibrium constant, B ₀, is not listed among the fixed parameters, because it varies somewhat withdifferent ionic conditions, different T, and different sequences. B ₀ for a given sequence under a particular condition can be obtained from the mean rise per bp or mean inverse bending elastic constant or mean inverse torsion elastic constant, the latter being the most sensitive and useful. The mean rise per bp is given by

(6)$$h = h_0- \delta /2 + f_b\delta $$

the mean inverse bending elastic constant is given by

(7)$$\left\langle {1/\kappa_\beta } \right\rangle = \lpar 1\,-\,f_b\rpar /\kappa _{\beta a} + f_b/\kappa _{\beta b}$$

and the mean inverse torsion elastic constant is given by

(8)$$\left\langle {1/\alpha } \right\rangle = \lpar 1\,-\,f_b\rpar /\alpha _a + f_b/\alpha _b$$

The acting torsion elastic constant is taken to be the reciprocal of the mean inverse torsion elastic constant and the acting bending elastic constant is the reciprocal of the mean inverse bending elastic constant. With knowledge of the model constants on the right-hand side (rhs) of Eqs. (6)–(8) and of the measured value on the left-hand side, f_b is readily calculated. For sufficiently long DNAs in the absence of tension, B ₀ is obtained from f_b by inverting Eq. (3) using the known constant value of J.

Variation of torsion elastic constant with T

Numerous FPA measurements on plasmid DNAs under different conditions in the absence of force (Delrow et al., Reference Delrow, Heath and Schurr1997; Delrow et al., Reference Delrow, Heath, Fujimoto and Schurr1998) are summarized as follows.

For ~100–110 mM univalent cations, P_f = 170 nm:

(9)$$ \alpha \lpar T \rpar \, = \,((6.1\, \pm \,0.4)- \lpar {0.0671} \rpar (T\,-\,293))\, \times \,10^{-19}J$$

For ~35–50 mM univalent cations, P_f = 180 nm:

(10)$$\eqalign{\alpha \lpar T \rpar = \lpar \lpar 6.2\,\pm \,0.4\rpar - \lpar {0.0671} \rpar \lpar T\,-\,293\rpar \rpar \,\times \,10^{{-}19}J}$$

For~39 mM univalent cations plus 5.5 mM Mg²⁺ ions, P_f = 150 nm:

(11)$$\alpha \lpar T \rpar = \lpar 8.9- \lpar {0.0614} \rpar \lpar T\,-\,273\rpar \rpar \,\times \,10^{{-}19}J$$

The ± values in Eqs. (9) and (10) are not standard deviations of measurements on a single DNA, but rather indicate the range of mean values for different DNAs. These values were corrected from measured lower bounds by using the relevant persistence lengths, P_f, reckoned for the fast-relaxing parts of the bending rigidities, as described elsewhere (Fujimoto and Schurr, Reference Fujimoto and Schurr1990; Schurr et al., Reference Schurr, Fujimoto, Wu, Song and Lakowicz1992; Okonogi et al., Reference Okonogi, Reese, Alley, Hopkins and Robinson1999; Naimushin et al., Reference Naimushin, Fujimoto and Schurr2000; Schurr, Reference Schurr2019a). Statistical uncertainties in these α(T) values are typically <~5% of the actual values.

Numerous experiments involving enzymes acting on DNA are conducted under ionic conditions similar to those for Eq. (11), but with 10 mM rather than 5.5 mM Mg²⁺. Because the addition of 5.5 Mg²⁺ to ~40 mM univalent cations increases α(T) by roughly 1.20-fold at 293 K, it might seem that adding 10 mM Mg²⁺ would raise α(T) a bit higher. However, the effect of added Mg²⁺ in the presence of sufficient univalent cations is believed to saturate at some relatively modest value, possibly 5.5 mM or less, so the enhancement due to 10 mM Mg²⁺ relative to 5.5 mM Mg²⁺ might not exceed 1.0 by very much. In this review, Eq. (11) is assumed to apply also for 10 mM Mg²⁺ in both 40 and 100 mM univalent cations.

It has been assumed that changes of ionic conditions among the three ranges in Eqs. (9)–(11) affect the properties of DNA primarily by altering f_b, while the properties of the individual a and b states remain invariant, and the same holds for variations in T in the 278–315 K range.

Intrinsic twists of the a and b states are nearly identical

The intrinsic twists (turns of one strand around the other per bp) of the a and b states must be practically the same in the low-force region under the prevailing conditions (~120 mM univalent cations, 296 K) for the following reasons. The observed linear increase in mean intrinsic twist with increasing force does not differ significantly from zero for tensions <1.3 pN (Gore et al., Reference Gore, Bryant, Nöllmann, Le, Cozzarelli and Bustamante2006). Even at 3.0 pN, the relative increase amounts to just one part in 2000. Any changes in mean intrinsic twist with increasing force from 0.13 to 0.92 pN are too small to produce observable changes in the reported zero-rotation positions under similar conditions (Mosconi et al., Reference Mosconi, Allemand, Bensimon and Croquette2009). The fraction f_b climbs from ~0.52 at 0.13 pN to ~0.91 at 1.3 pN, but that evidently causes no detectable change in intrinsic twist.

Intrinsic twists of the a and b states decline similarly with increasing T

Measurements of the intrinsic twist, t ₀ (turn bp⁻¹ K⁻¹), of single DNAs held at various tensions (0.3, 0.5, 0.7, 0.8, and 0.9 pN) at temperatures in the range of 296–315 K decreased linearly with slopes, dt ₀/dT = − 3.19 × 10⁻⁵ and− 2.92 × 10⁻⁵ turn bp⁻¹ K⁻¹ for, respectively, 20.6 and 7.9 kbp constructs in 160 mM univalent cations (Kriegel et al., 2018). These slopes were largely independent of the applied tension, although a marginally significant upward trend appeared at 0.8 and 0.9 pN.

With increasing T from 273 to 299 K, the intrinsic twist of large circular DNAs decreased linearly with slopes, dt ₀/dT = − 3.389 × 10⁻⁵ turn bp⁻¹ K⁻¹, for 5300–9850 bp DNAs, and, dt ₀/dT = − 3.14 × 10⁻⁵ turn bp⁻¹ K⁻¹, for a 2200 bp DNA under the same ionic conditions (2 mM Mg²⁺, 10 mM Tris HCl, pH 8.0, 1 mM Na₃EDTA) (Depew and Wang, Reference Depew and Wang1975). Much smaller ~200 bp circular DNAs exhibited similar slopes, dt ₀/dT = − 3.19 × 10⁻⁵ turn bp⁻¹ K⁻¹, over the 278–315 K range under slightly different ionic conditions (~40 mM univalent cations plus 10 mM Mg²⁺) (Geggier et al., Reference Geggier, Kotlyar and Vologodskii2011). Because such small circles possess significant tension and a correspondingly higher fraction of the bp in the b state, these results suggest that the intrinsic twists of the a and b states are not only similar at ambient T, but also vary with T in a similar manner.

Mean rise per bp and persistence length at zero force

In zero force, this value of α ₀ yields, f_{b 0} = 0.464, B ₀ = 0.999357, and h = 0.323 nm, which coincides with the SAXS results on 10–35 bp DNAs labeled with gold colloids under the same conditions (Zettl et al., Reference Zettl, Mathew, Seifert, Doniach, Harbury and Lipfert2016). It also yields, P ₀ = h $\kappa _\beta$/k _BT = 52.0 nm, which lies in the canonical range, 50 ± 4 nm.

Relative extension versus force

The following example demonstrates how to use the model in a calculation. One begins with the known force, F, and Eq. (3) with B_F in place of B. After setting, 〈cos θ〉 = z_F/L_F in Eq. (4), Eq. (5) is used to express z_F/L_F in terms of F, T, and P_F = h_F κ_βF/k _BT. Both h_F and κ_βF are expressed in terms of f_b and constants of the model using Eqs. (6) and (7). Now B_F is expressed entirely in terms of F, T, constants of the model, and f_b, so substitution of Eq. (4) into Eq. (3) yields a transcendental equation with f_b as the only unknown. Mathematica is used to solve numerically for f_b, from which the values of h_F, κ_βF, P_F, and z_F/L_F are readily evaluated for each value of the force. The experimental data are in the form, z_F/z _max, where the value of z _max is unknown. Equation (5) is re-written as

(12)$$\;z_F/z_{\max } = \lpar L_{\max }/z_{\max }\rpar \lpar L_F/L_{\max }\rpar \lpar 1\,-\,\lpar 0.5\rpar /\lpar P_FF/k_{\rm B}T\,-\,1/32\rpar ^{1/2}\rpar $$

The quantity, L _max = N(0.325 + 0.0487/2) = N(0.349) nm is the largest possible value of the contour length L_F attainable by the two-state model_. The value, z _max = 0.975 L _max, was adopted, so that the theoretical value, z_3.9/z _max = 0.947, precisely matches the corresponding experimental value for F = 3.9 pN (Mosconi et al., Reference Mosconi, Allemand, Bensimon and Croquette2009). This seemingly arbitrary protocol, which improves the fit of the model to the experimental data, is actually self-consistent, as will be shown. In any case, it is expected that z _max will be slightly less than L _max at the highest force in the data set, because z _max was obtained by fitting a non-stretchable WLC model to data over the whole range of forces, including lower forces, where the DNA has somewhat shorter contour lengths. The fraction of the bp in the longer b state at F = 3.9 pN is f_b = 0.989, so the DNA is practically stretched to its maximum rise per bp, h _max = 0.349 nm, at that point, even though it is not quite fully extended (by orientation). The experimental and theoretical values of z_F/z _max are compared in Fig. 2. The agreement is rather good. Relative deviations from the experiment by as much as 1.4% and 1.0% occur only at the two lowest forces (0.18 and 0.25 pN), and likely lie within the experimental errors of the measured values.

Fig. 2. Relative extension, z_F/z _max, versus pulling force F for an untwisted DNA at 293 K. Black disks are theoretical values predicted by the two-state model with the parameters in Table 1, and gray disks are experimental values of Mosconi et al. (Reference Mosconi, Allemand, Bensimon and Croquette2009). The theoretical values were computed for the two-state model by assuming that, z _max = (0.975) L _max. Where gray experimental points are not visible, they lie underneath black theoretical points. The dashed curves were computed for non-stretchable WLC models with different values of P and the same z _max. From top to bottom, the dashed curves correspond to P = 90, 50, 44, and 38 nm. The curve for 44 nm has the smallest chi-squared value of any of the dashed curves.

The assumed value, z _max/L _max = 0.975, is self-consistent in the following sense. The mean extension per bp, z ₁ = z_F/N = z_F h_F/L_F = h_F z_F/L_F, of the two-state model in the z-direction was reckoned at each experimental value of the force, and these were regarded as ‘measured’ values. The non-stretchable WLC model with P and h as adjustable parameters was fitted to these ‘measured’ z ₁-values. A minimum in the sum of squared differences between the predictions of the non-stretchable WLC model and the ‘measured’ values (from the two-state model) was found at P = 47 nm and h = 0.341 nm. The latter value is the maximum extension per bp, z1 _max, of the best-fit non-stretchable WLC model. The ratio of the maximum extension of the non-stretchable model to that of the two-state model is z1 _max/h _max = 0.341/0.349 = 0.977, which is practically identical to the value 0.975 assumed above. Hence, the two-state model accurately mimics the experimental data, and the assumed value, z1 _max/h _max = 0.975 is self-consistent.

A second minimum in the sum of squared deviations between predictions of the non-stretchable WLC model and the predictions of the two-state model was found at P = 38 nm and h = 0.351 nm. This minimum was rejected because the values obtained must apply at zero force, as well as at forces in the experimental range, 0.18–3.9 pN. The value, P = 38 nm, is much too small for an unstressed DNA in 0.1 M univalent cations at 293 K, and the value, h = 0.351 nm, significantly exceeds the mean rise per bp, 0.323 ± 0.01 nm, obtained from the X-ray scattering data.

The dashed curves in Fig. 2 indicate how poorly the non-stretchable WLC models with various choices of persistence length fit the experimental relative extension values.

How cooperativity works

Set B ₀ = 1.0 and x = (Fδ〈cos θ〉/k _BT). When x ² << 4 J ² << 1.0, Eq. (3) can be written to first order in x as follows:

(13)$$f_b = \lpar {1/2} \rpar \lpar 1 + x/\lpar {2J} \rpar \rpar $$

Using δ = 0.0487 nm, F = 0.18 pN, and the experimental value 〈cos θ〉 = z_F/L_F = 0.63, the values x = 1.36 × 10⁻³, x ² = 1.86 × 10⁻⁶, J = 4.4855 × 10⁻³, and 4J ² = 8.05 × 10⁻⁵ are obtained. Hence, 4J ² exceeds x ² by a factor of 42, and the conditions for validity of Eq. (13) hold. Since 1/2J = 111, the value of x has been effectively enhanced by 111-fold due to the smallness of the cooperativity parameter. In effect, the difference in length for 111 bp (i.e. 111 × δ), instead of that for a single bp, governed the response of the equilibrium to the force. The fraction of the bp in the b state has jumped from 0.50 in the absence of force to 0.58 in the presence of the small force, despite the fact that x = Fδ〈cos θ〉/k _BT = 1.36 × 10⁻³ is so small.

A key aspect of the two-state model is that the properties of the individual a and b states are assumed to remain practically invariant to small perturbations of any kind, as would be expected for a ‘solid’ material. The observed variations are attributed entirely to shifts in the a⇔b equilibrium. The surprisingly great sensitivity of DNA properties to modest changes in ionic conditions, temperature, pulling, twisting, bending, and sequence changes is attributable to two aspects of the model: (1) the very great cooperativity arising from the very small value of J and (2) the fact that B ₀ values are so close to 1.000, such that (B ₀ − 1)² << 4B ₀J ². Only in such a circumstance is it possible for the large cooperativity to maximally amplify the effect of small perturbations.

Effective torsional rigidity and intrinsic torsion elastic constant versus force

The effective torsional rigidity, C _eff, of a DNA under tension is less than its intrinsic value because the helix-axis absorbs some of the imposed twists into writhe at a lower energy cost than would be required to absorb it all as a pure twist. Dividing the measured torque by the imposed twist angle per unit length yields the experimental value of C _eff. The theoretical value of C _eff is reckoned using the model parameters in Table 1 and α ₀ = 6.3 × 10⁻¹⁹ J, which differs only slightly from the value, α ₀ = 6.1 × 10⁻¹⁹ J, adopted previously (Schurr, Reference Schurr2019a, Reference Schurr2019b). MN torque theory expresses the reduced (by k _BT) torque in terms of the tension, F; the (small) excess twist angle per unit length; the intrinsic torsional rigidity, C _in = h_F α_F; the persistence length, P_F; and the reduced torque. Explicit expressions are presented elsewhere (Moroz and Nelson, Reference Moroz and Nelson1998; Schurr, Reference Schurr2015). Such an expression was solved for the reduced torque for every experimental value of the force by using Mathematica, and each value of C _eff was obtained from that. As part of the same calculation, the values of f_bF, h_F, P_F, α_F, and the intrinsic torsional rigidity, C _inF = α_Fh_F, were also reckoned for each force. The results are shown in Fig. 3a. The experimental C _eff values fall ever further below the predictions of the non-stretchable TWLC model with decreasing force, demonstrating the inadequacy of any model with a torsion elastic constant that is independent of the tension. The two-state model provides a reasonable compromise fit to three different data sets extending from 0.20 to 6.5 pN. Moreover, it yields, C _eff = 411 × 10⁻³⁰ Jm, for a force of 15 pN, in agreement with the experimental values, C _eff = (410 ± 30) × 10⁻³⁰ and (440 ± 40) × 10⁻³⁰ Jm at the same force (Bryant et al., Reference Bryant, Stone, Gore, Smith, Cozzarelli and Bustamante2003).

Fig. 3. (a) Effective torsional rigidity, C _eff, versus force, F, for a single slightly twisted DNA molecule held at force, F. Thick black curve represents the theoretical values computed using the two-state model parameters given in Table 1 at T = 293 K, and α ₀ = 6.3 × 10⁻¹⁹ J. The gray squares are experimental C _eff data of Mosconi et al. (2009), and the black disks are from two different magnetic torque tweezers experiments of Lipfert et al. (Reference Lipfert, Kerssemakers, Jager and Dekker2010) and Lipfert et al. (Reference Lipfert, Wiggin, Kersemakers, Pedaci and Dekker2011). The dashed black curve is computed for a non-stretchable TWLC with P = 50 nm, a rise per bp, h = 0.34 nm, and an intrinsic torsion elastic constant, α = 12.32 × 10⁻¹⁹ Jm. (b) Predicted fraction, f_b, of bp in the b state (left axis, thick black curve) and α/α _b (right axis, dashed curve) versus force, F, for the same DNA as in (a). This curve was computed for the same two-state model with the same parameters as in (a).

The computed fraction, f_b, of the bp in the b state and (left axis) and the ratio, α/α_b (right axis) are plotted versus F in Fig. 3b. Straightforward extrapolation to F = 0 indicates that neither f_b nor α/α_b declines below about 0.5 in the absence of force for this particular DNA.

An alternative single-state model with a potential energy function that included twist–bend coupling and anisotropy of bending terms, in addition to pure bending and twisting terms, was proposed to account for the experimental data (Nomidis et al., Reference Nomidis, Kriegel, Vanderlinden, Lipfert and Carlon2017). However, in subsequent work that suggestion was effectively withdrawn, because it was based upon simulations using a too long elementary subunit length. A deeper analysis revealed that regardless of the local potential function, the behavior over long distances (several turns or more) becomes like that of a TWLC with an effective torsional rigidity, but no twist–bend coupling or anisotropy of the bending. This is due to the rotational averaging of the properties over several turns. The long-distance torsional rigidity depends upon the local twist–bend coupling and bending anisotropy potential parameters in the derived manner (Nomidis et al., Reference Nomidis, Skoruppa, Carlon and Marko2019). However, this new theory is evidently incapable of fitting the experimental data at both low and high forces with a single set of parameters. That leaves the present two-state cooperative model as the only currently viable option for rationalizing these single-molecule torque data.

A comparison of linear and circular forms of a 181 bp DNA

The observed ~1.5-fold enhancement of the torsion elastic constants of two different samples of a 181 bp DNA upon circularization established conclusively that a coherent bend of 2°/bp significantly elevated the torsion elastic constant (Heath et al., Reference Heath, Clendenning, Fujimoto and Schurr1996). The coherent bending strain seemed to be the most likely cause of the increased torsion elastic constants in the circular DNAs. However, even higher torsion elastic constants than those in the 181 bp circles prevailed in DNAs held at 15 pN tension, where bending was almost entirely suppressed (Bryant et al., Reference Bryant, Stone, Gore, Smith, Cozzarelli and Bustamante2003). Was the torsion elastic constant actually a constant, as many people believed, or was the tension causing the DNA to raise its intrinsic torsion elastic constant from a lower to a higher value already in the low-force regime? That possibility was affirmed by applying the MN theory to experiments in that low-force regime, as demonstrated above. Increasing tension from 0.2 to 4.0 pN increased the intrinsic torsion elastic constant from values comparable to those of unstrained linear DNAs toward saturation near 3.5–4.0 pN at a plateau value close to the previously measured values at 15 pN (Schurr, Reference Schurr2015).

A year and a half later it finally dawned on me that bending a rod with a non-vanishing bending rigidity into a circle necessarily introduces tension into the rod. The basic argument is simple. The strain energy, U, of a thin elastic rod bent into a planar circle is proportional to (i) the length of the rod, L, (ii) the bending rigidity, A, and (iii) the square of the curvature. The curvature is proportional to 1/R, where R = L/(2π) is the radius of the circle, so U is ultimately proportional to 1/L. Because U decreases with increasing L, it acts like a tension to stretch the length of the DNA. More precisely, the axial tension acting to stretch the DNA is given by

(14)$$ F_{{\rm str}}\, = \,\,-\,{\rm d}U/{\rm d}L_{{\rm cir}}\, = \,(A_{{\rm cir}}/2) \lpar {2\pi } \rpar ^2/L_{{\rm cir}}^2$$

where A _cir = κ_{β cir}h _cir, L _cir = Nh _cir, and h _cir is the rise per bp, N is the number of bp, and the subscript ‘cir’ applies to the circular form of the DNA.

Subsequent to the original paper on this topic (Schurr, Reference Schurr2017), I discovered that the number of bps, N = 181, is not large enough to render the λ ₋^N term from the lower eigenvalue entirely negligible in comparison with the λ ₊^N term in Eq. (1). Consequently, that term should be included in the expression for χ. Using the symbolic derivative in 〈n_b〉 yields a considerably more complex algebraic expression for the rhs of Eq. (3) in terms of B_F and J. Apart from that, the calculation of f_b is still straightforward for the linear 181 bp DNA, beginning with α ₀ = 6.315 × 10⁻¹⁹ J for the corresponding long (N > 1000 bp) linear DNA. The value, B ₀ = 0.999357, applicable to both long and short unstrained DNAs, is obtained by inverting Eq. (3) using the simple expression valid for large N.

The situation is somewhat more complex for the circular 181 bp DNAs. For a circular DNA, the last bp must be in the same state as the bp that precedes the first bp, in which case the partition function becomes, Tr|M^N| = λ ₊^N + λ ₋^N. The value, B ₀ = 0.999357, is the same as for the linear DNA. B_F can be expressed as in Eq. (4), but with two modifications. First, 〈cos θ〉 = 1.0, since the axial stretching force is always directed along the local helix-axis. Second, the linear DNA contains 17.238 turns, whereas the predominant topoisomer contains precisely 17 turns, and is underwound by 0.238 turns or 1.496 radians. The energy per bp arising from this incorporated twist (−1.496/181 rad/bp) is expressed in terms of α_a, α_b, T, and the acting value of α (Schurr, Reference Schurr2017; Supporting Information), and α, in turn, is expressed in terms of f_b, α_a, and α_b via Eq. (8). This twisting energy appears as an additional exponential factor analogous to the exponential stretching factor in Eq. (3), but with the opposite sign in the exponent. It has a comparatively minor effect. The stretching force, F _str, is expressed in terms of N, T, and the unknown κ_{β cir}, and h _cir, both of which are expressed in terms of model constants and f_b. Now, B_F is expressed entirely in terms of model constants, T, B ₀, and f_b. Upon substituting B_F into the modified, now more complex, rhs of Eq. (3), it becomes a transcendental equation for f_b that is solved numerically by Mathematica. With f_b in hand, all properties are readily calculated.

The theoretical results for the 181 bp linear and circular DNAs, and for the corresponding long linear DNA are compared with each other and with experimental values in Table 2. Agreement of the two-state model with experiment is remarkably good for both linear and circular 181 bp DNAs. None of these experimental values were used in the original fitting protocol to determine the model constants. A comparison of theoretical values for the linear 181 and long linear DNAs indicates that inclusion of the contribution of the second eigenvalue has only a minor effect in the case of the linear 181 bp DNA.

Table 2. Comparison of theoretical values for 181 bp linear, 181 bp circular, and long linear DNAs in 0.1 M univalent cations at 293 K with each other and with FPA experiments on 181 bp DNAs of the same sequence (Heath et al., 1996)

^a Torsion elastic constant of the effective springs between bps reckoned from the average value of 1/α.

^b Bending elastic constant of the effective springs between bps reckoned from the average value of 1/κβ.

^c Average rise per bp.

^d Persistence length reckoned from P and hav.

^e Average fraction of the bp in the b state.

Supercoiling free energies for ~210 bp circular DNAs

The free energy difference between two topoisomers of a circular DNA with N  bp is believed to be expressible as

(15)$$\Delta A = A\lpar n \rpar - A\lpar p \rpar = k_{\rm B}T\lpar E_{\rm T}/N\rpar \lpar \Delta n^2- \Delta p^2\rpar $$

where E _T is the twist energy parameter, Δn = n − t ₀ and Δp = p – t ₀ are the linking differences in turns, n and p are the invariant integral linking numbers (turns of one strand around the other) in the circles, and t ₀ is the total intrinsic twist (turns), which is generally non-integral. This relation was found empirically (Depew and Wang, Reference Depew and Wang1975; Pulleyblank et al., Reference Pulleyblank, Shure, Tang, Vinograd and Vosberg1975) for long DNAs, and proved analytically for a discrete TWLC model, whenever the distribution of writhe in the corresponding unconstrained DNA is Gaussian (Wu et al., 1988). (An unconstrained DNA does not obey the Calugareanu, White, Fuller topological constraint, according to which, n = t + w, must be an integer, where t (turns) is the normally non-integral total twist (turns), and w is the normally non-integral writhe (turns). Its writhe is reckoned by taking no account of its twisting coordinates, and its net twist, t – t ₀, is reckoned by taking no account of its bending coordinates, or the path of the helix-axis.) A simple rearrangement of the expression obtained for E _T by Wu et al. (Reference Wu, Song, Clendenning, Fujimoto, Benight and Schurr1988) and conversion from units of radians to turns yields

(16)$$E_{\rm T} = \lpar {N/2} \rpar /\lpar \left\langle {w^2} \right\rangle _u + \left\langle {{\lpar t- t_0\rpar }^2} \right\rangle _u\rpar $$

where 〈(t – t ₀)²〉_u = (1/(2π)²)Nk _BT/α is the variance of the net twist (turns²) of the unconstrained DNA, 〈w ²〉_u = (1/(2π)²)Nk _BT/κ_w is the variance of the Gaussian distribution of writhe (turns²) for the unconstrained DNA, and κ_w is the unknown effective Hooke's law constant for the writhe of the unconstrained DNA. The variance, 〈(t – t ₀)²〉_u, is calculated directly from the acting α according to the expression above. Although κ_w is unknown, the variance of the writhe was simulated for a wide range of lengths, persistence lengths, and effective helix diameters, and the results were combined to yield comparatively simple and useful analytical expressions (Geggier et al., Reference Geggier, Kotlyar and Vologodskii2011) that can be used to reckon〈w ²〉_u from N, h, κ_β, and T plus the effective helix diameter, d, which is intended to account for self-repulsions due to electrostatic interactions. Equations (15) and (16) are valid, when, and only when, the distribution of writhe in the unconstrained DNA is Gaussian. That circumstance is expected to prevail for sufficiently small circles with linking differences well below the midpoint (~1.3) of the first transition in writhe, and again for large circular DNAs at any thermally accessible value of the linking difference. In the intermediate region, the distribution of writhe is not expected to be Gaussian (Gebe and Schurr, Reference Gebe and Schurr1996).

Acting parameters of the circle were reckoned as done above for the 181 bp circular DNA, except that α ₀ was taken from Eq. (11) instead of (9), because of the different prevailing conditions. E _T was reckoned from 〈w ²〉_u and 〈(t − t ₀)²〉_u via Eq. (16). Values of 〈w ²〉_u were calculated using Eqs. (2) and (3) of Geggier et al. (Reference Geggier, Kotlyar and Vologodskii2011) with d = 6 nm (Schurr, Reference Schurr2019a). The results for each length of DNA are presented in Table 3.

Table 3. Calculated properties of circular DNAs with 205–4362 bp in ~40 mM univalent cations plus 10 mM Mg2+ at 310 K.

Beginning with α _lin for a long linear DNA reckoned from Eq. (11), values were computed as described above for the 181 bp DNA. E _T is reckoned from Eq. (16) and compared with experimental values obtained from topoisomer distributions (Horowitz and Wang, Reference Horowitz and Wang1984) (HW). The dual values of the experimental E _T and helix repeat, H, were obtained by least-squares fitting of the experimental TRs for 205, 207, and 217 bp DNAs to Eqs. (6b) and (5), respectively, of HW.

For any small ~200 bp circular DNA, at most two topoisomers have sufficiently similar populations that their ratio can be reliably assessed, and that requires that its linking difference be not too far from half-integral (where the two topoisomers would have the same free energy and be equally abundant). The ratio of the topoisomer population with the greater linking number to that with the smaller linking number was measured for each of three circular DNAs containing 205, 207, and 217 bp of largely common sequence in ~30 mM univalent cations plus 10 mM Mg²⁺ at 310 K (Horowitz and Wang, Reference Horowitz and Wang1984) (HW). The negative logarithm of such ratios yields the difference in free energy between the topoisomer with the larger integral linking number and that with the next smaller value, divided by k _BT. Least-squares fitting of Eq. (5) of HW to the measured TRs for each of the three DNAs simultaneously yielded, E _T = 3980, and helix repeat, H = 10.54 bp/turn, but the fitting of Eq. (6b) of HW to the ratios of the TRs for two different DNAs for all three pairs of DNAs simultaneously yielded, E _T = 3910, and helix repeat, H = 10.46. This difference in results of fitting one formula for individual ratios and fitting another for the ratios of individual ratios for pairs of DNAs is due to the fact that neither of the fits with the global best-fit parameters provides satisfactory fits of the individual TRs of all three DNAs. HW acknowledged the sizeable fitting discrepancies for the individual DNAs by admitting a range of E _T from 3700 to 4200 for a 210 bp DNA. In any event, the value, E _T = 3867, predicted by the two-state model for a 210 bp DNA underestimates the reported, E _T = 3910 and 3980, by about 1% and 3%, respectively, which is well inside the experimental uncertainty.

The theoretical and experimental E _T values for 551 bp DNAs disagree substantially, as expected because neither Eq. (15) nor (16) is valid in this intermediate region, where the distribution of writhe in the unconstrained DNAs is not Gaussian. Equation (15) for a simulated 468 bp circular DNA becomes invalid for linking differences >0.5 turn (Gebe and Schurr, Reference Gebe and Schurr1996). This implies a non-Gaussian distribution of writhe in the unconstrained DNA, which arises from thermally accessible Fig. 8 configurations in addition to slightly writhed circles. The theoretical E _T value is somewhat closer to the experimental value for the 1188 bp DNA and exceeds by only 1.02-fold, E _T = 1220, for the 2434 bp DNA, and by 1.01-fold, E _T = 1130, for a 4362 bp (pBR322) DNA. All of these E _T measurements have ~6–10% relative statistical errors. Another experimental value, E _T = 1155 ± 65, was also reported for a 4363 bp pBR322 DNA under similar conditions (Naimushin et al., Reference Naimushin, Clendenning, Kim, Song, Fujimoto, Stewart and Schurr1994), so the theoretical value actually lies between the two reported experimental values. In the range of circle sizes where it is expected to apply, the two-state cooperative model provides an accurate account of the E _T values.

Variation of the bending rigidity with T

A recent review of the relevant literature concluded that the effective bending rigidity of DNA, reckoned from its persistence length, actually undergoes a modest net increase with increasing T over the range from 278 to 315 K (Schurr, Reference Schurr2019a). This conclusion was derived from experimental results pertaining to long linear bacteriophage and plasmid DNAs, linear restriction fragments thereof, synthetic linear DNAs of various lengths, and a relaxed 2686 bp circular DNA. The experimental results were obtained by a variety of methods: sedimentation velocity (Gray and Hearst, Reference Gray and Hearst1976), dynamic light scattering (DLS) (Wilcoxon and Schurr, Reference Wilcoxon and Schurr1983), off-field decays of the electric dichroism (Diekmann et al., Reference Diekmann, Hillen, Morgeneyer, Wells and Pörschke1982; Pörschke, Reference Pörschke1991), and electric birefringence (Lu et al., 2001, Reference Lu, Weers and Stellwagen2002), and two different tethered particle distribution measurements (Driessen et al., Reference Driessen, Sitters, Laurens, Moolenaar, Wuite, Goosen and Dame2014; Brunet et al., Reference Brunet, Salomé, Rousseau, Destainville, Manghi and Tardin2017). The latter study revealed a time-average blurring artifact of the method, for which the authors corrected both their own data and those from the earlier study. After such corrections, all of these studies indicated that the bending rigidity underwent a net increase or, in the case of one sample, did not change, as T was varied from 278 to 315 K. The opposite conclusion was obtained in the original analysis of the measured topoisomer distributions and E _T values for a 2686 bp circular pUC19 DNA (Geggier et al., Reference Geggier, Kotlyar and Vologodskii2011). However, that analysis was based upon the invalid assumption that the torsion elastic constant is independent of temperature and exhibits the value, α = 9.12 × 10⁻¹⁹ J, under all conditions. When the same E _T versus T data were analyzed using the measured T-dependent values of α in Eq. (11), the bending elastic constants obtained also exhibited a modest increase over the range from 278 to 315 K, as shown in Table 4. Particularly noteworthy are the opposite directions of the changes in the torsion and bending elastic constants with increasing T, which are unexpected and in need of explanation.

Table 4. Pertinent quantities for pUC19 DNAs (N = 2686 bp) in buffers containing 39.9 or 37.2 mM univalent cations plus 10 mM Mg2+ at various T. The E_T values were obtained from topoisomer distributions by Geggier et al. (Reference Geggier, Kotlyar and Vologodskii2011) and 〈Lk ²〉_u ≡ (〈w ²〉_u + 〈(t – t0)²〉_u) from N/(2E _T)

The α values come from Eq. (11), and the f_b and h are reckoned from Eqs. (8) and (6), respectively. The P values were determined by trial and error using the protocol of Geggier et al. (Reference Geggier, Kotlyar and Vologodskii2011), but with the length computed from N and h. Then, κ_β = k _BTP/h.

Several of the cited experimental studies revealed a significant sequence dependence of both the magnitude of the effective bending rigidity and its variation with temperature. The effective bending rigidity, reckoned from the persistence length, is a complex quantity with contributions to the mean-squared bending angle, or inverse persistence length, from rapidly relaxing bends, non-relaxing (permanent) bends, and slowly relaxing bends (Naimushin et al., Reference Naimushin, Fujimoto and Schurr2000).

Fitting Eq. (7) to the κ_β values in Table 4 by using the f_b values in that same table yielded the best-fit values of κ_βa and κ_βb, which are the model parameters in Table 1. The fit of 1/κ_β to f_a = 1 − f_b is reasonably good, as demonstrated in Fig. 4 (top axis). The question arises, whether the two-state model could fit any set of κ_β values that varied smoothly with T. Several tests were conducted with assumed linear, quadratic, or hyperbolic variations of κ_β with increasing T, but only relative variations close to that in Table 4 gave fits as good as that in Fig. 4. From that perspective, this good fit is a fairly strong argument in favor of a two-state model, wherein the change in elastic properties is attributed to a shift in the relative occupations of the two states, while the properties of each state remain invariant to changes in T over the limited range from 278 to 315 K.

Fig. 4. Inverse bending elastic constant (1/κ_β) versus fraction of bp in the a state (f_a) (top axis), and versus inverse torsion elastic (1/α) (bottom axis). The black circles represent the values in Table 4, which were derived by Schurr (Reference Schurr2019a) from experimental E_T values of Geggier et al. (Reference Geggier, Kotlyar and Vologodskii2011).

A strong linear anti-correlation between the inverse bending and inverse torsion elastic constants is also demonstrated in Fig. 4 (bottom axis). This argues that changes in both elastic constants stem from a common cause, such as the T-dependent shift in the relative populations of two states with different values of the torsion and bending elastic constants. With increasing T the shift is toward the shorter, torsionally softer, flexurally stiffer a state.

Is the two-state cooperative transition the ‘broad pre-melting transition’?

Beginning in the 1960s, several groups reported significant changes in the UV absorbance, UV CD, UV optical rotation, infrared absorption, and Raman vibrational spectra of many, but not all, DNAs, as T was increased from 0 °C up to the melting region (Brahms and Mommaerts, Reference Brahms and Mommaerts1964; Luck et al., Reference Luck, Zimmer, Snatzke and Söndergrath1970; Gennis and Cantor, Reference Gennis and Cantor1972; Palecek and Frick, Reference Palecek and Frick1972; Ramm et al., Reference Ramm, Vorobev, Birshtein, Bolotina and Volkenstein1972; Shie et al., Reference Shie, Kharitonenkov, Tikhonenko and Chirgadze1972; Studdert et al., Reference Studdert, Patroni and Davis1972; Sarocchi and Guschlbauer, Reference Sarocchi and Guschlbauer1973; Usatyi and Shlyakhtenko, Reference Usatyi and Shlyakhtenko1973; Rimai et al., Reference Rimai, Maher, Gill, Salmeen and McCormick1974). The amplitude of the positive CD band (ellipticity per mole of bases) near 273 nm increased with increasing T, whereas that of the negative band near 247 nm remained more or less unchanged. In contrast, in the melting region, the amplitudes of both bands decreased abruptly for all DNAs, including those that showed negligibly small amplitudes of pre-melting change. Observed changes below the melting region were ascribed to a ‘broad pre-melting transition.’

T7 DNA exhibited a significant amplitude of pre-melting change in its CD spectra, and also in its UV difference absorption spectra, as noted earlier (cf. Fig. 1). With increasing T, the absorbance of T7 DNA between 239 and 262 nm declined, but beyond 262 nm it rose. (The absolute absorption spectrum reaches a maximum near 262 nm.) If the two states monitored correspond to a and b, then these observations imply that the absorbance band of the a state, favored at higher T, is slightly red-shifted from that of the b state, favored at lower T. This suggests that unexcited bases in the a state are in some sense farther apart and/or more weakly interacting than in the b state, and thus better able to accommodate the expanded electron cloud of the excited base at a lower energy cost than is the case for the b state. Greater separations and/or weaker interactions with neighbors would be consistent with the greater entropy of the a state, which is responsible for its increasing prevalence with increasing T (Schurr, Reference Schurr2017; Supporting Information II).

Interest in pre-melting changes waned in the late 1970s, in part because evidence of significant changes in other physical properties of DNA was lacking.

A much later study found that pre-melting transitions in the CD spectra of three different DNAs, including an 1876 bp HaeII restriction fragment of pBR322, were accompanied by substantial declines in torsion elastic constant, α, with increasing T (Delrow et al., Reference Delrow, Heath, Fujimoto and Schurr1998). Moreover, DLS studies of that 1876 bp DNA at large scattering vectors indicated that its dynamic bending elastic constant, κ_β, increased significantly from 293 to 333 K in the pre-melting range. These findings suggest that increasing T favors the a state in these DNAs. Addition of 0.5 M tetramethylammonium chloride, the cation of which is known to bind preferentially to AT-rich regions, substantially diminished these pre-melting changes in CD, α, and κ_β with increasing T. This implies that the cation, TMA⁺, prefers to bind to AT-rich regions in the b state, rather than to those in the a state, thereby stabilizing those regions in the b state, and inhibiting their transition, and that of their neighbors, toward the a state with increasing T.

At first glance, the wide T-range of the pre-melting transition seems to contradict the notion of a highly cooperative transition. However, according to Table 4, the fraction of the bp in the b state reckoned according to the model declines steadily from 0.755 to 0.465, as T is increased from 278 to 315 K. Over the same T range, the value of B reckoned by inverting Eq. (3) declines from 1.0053 to 0.999371, and its associated standard state-free energy difference (in units of k _BT) between the b and a states of a single bp, ΔG ⁰/k _BT = − ln(B), declines from +5.31 × 10⁻³ to −6.29 × 10⁻⁴. Such exceedingly small changes are manifested only because of the very great cooperativity of the transition. In such a circumstance, the transition can take place gradually over a wide T range.

If the pre-melting increase in CD at 273 nm reflects an increase in the fraction of bp in the a state with rising T, consistent with Fig. 4 for the two-state model and Eq. (9), then a decrease in CD at 273 nm should accompany any perturbation that shifts the a⇔b equilibrium toward the b state. Circularizing the aforementioned 181 bp DNA, which increased its α by ~1.5-fold (presumably by shifting the equilibrium toward the b state) also decreased its CD at 273 nm by 0.9-fold (Heath et al., Reference Heath, Clendenning, Fujimoto and Schurr1996) as expected.

Range of validity of the two-state cooperative model

The range of temperatures considered for the two-state model (278–315 K) is considerably smaller than the full pre-melting range (273–333 K). The limited range adopted here is not fundamental, but instead was invoked to ensure approximate constancy of the model parameters, such as the cooperativity parameter J, the rises per bp, h_a = h ₀ – δ/2 and h _b = h ₀ + δ/2, and the elastic constants, α_a, α_b, κ_βa, and κ_βb, of the individual states, over the T range considered. A two-state equilibrium may well prevail over a larger T range, but J and properties of the individual states are likely to vary significantly over a wider T range.

Why two states of very similar free energy instead of just one?

Duplex DNA can be regarded as a four ‘strand’ helical polymer with two different central columns of stacked complementary bases interconnected by their H-bonds, plus two outer covalently connected sugar-phosphate backbones. The two central interconnected columns likely experience some frustration of their own. If the backbones, which are very rigid in some respects, are sufficiently incommensurate with the attachment configurations on the two base columns, then significant overall frustration would likely prevail in any final structure. In such a case, one might expect to observe at equilibrium two or more distinct structures of similar energy, which are also incompatible in the sense that the free energy required for them to be nearest neighbors is extremely large. Such a frustrated material with large domains would be expected to display slow relaxations of certain strains following changes in their corresponding stresses. Indeed, slow relaxation of the UV difference spectra upon changing the sample temperature over the range 3–20 °C was already mentioned in the discussion of Fig. 1. Numerous examples of slow (hours, days, or weeks) changes in the torsion elastic constant, CD spectrum, DLS at large scattering vector, and even electrophoretic mobilities subsequent to changes in superhelical stress and other environmental factors were reported by our lab (Shibata et al., Reference Shibata, Wilcoxon, Schurr and Knauf1984; Langowski et al., Reference Langowski, Benight, Fujimoto, Schurr and Schomburg1985; Wu et al., Reference Wu, Song, Clendenning, Fujimoto, Benight and Schurr1988; Wu and Schurr, Reference Wu and Schurr1989; Song et al., Reference Song, Fujimoto, Wu, Thomas, Shibata and Schurr1990; Wu et al., Reference Wu, Fujimoto and Schurr1991; Schurr et al., Reference Schurr, Fujimoto, Wu, Song and Lakowicz1992; Kim et al., Reference Kim, Fujimoto, Furlong, Sundstrom, Humbert, Teller and Schurr1993; Naimushin et al., Reference Naimushin, Clendenning, Kim, Song, Fujimoto, Stewart and Schurr1994; Schurr et al., 1998; Brewood et al., Reference Brewood, Delrow and Schurr2010).

Sequence preferences of the a and b states of the two-state model

By assuming that the pre-melting transition corresponds to a shift of the a⇔b equilibrium toward the a state with increasing T, information about the preferences of different types of sequences for a or b states might be gained. A larger fraction, f _G+C, of G + C bases in the DNA was modestly correlated with smaller increases in the difference spectrum at 287 and 295 nm for a given rise in T (SG). Hence, a larger f _G+C favors a smaller f_b and a greater f_a at any given T, because only those bp in the b state can convert to the a state with rising T, thereby increasing the difference spectrum. These same changes in the difference spectrum at both 287 and 295 nm were found to correlate better with quotients, q _like = (f _TpT + f _CpC)/f _G + C and q _all = (f _TpT + f _CpC + f _TpC + f _CpT)/(2f _G + C), where the fractions in the numerator pertained to ‘like’ (TpT and CpC) and ‘all’ (TpT, CpC, TpC, and CpT) pyrimidine–pyrimidine dinucleotide steps out of all possible steps in the DNA. Such nearest-neighbor frequencies were the only sequence information available at that time. The convention, y = pyrimidine, r = purine, is adopted here. Evidently, increasing the fraction of yy dinucleotide steps, while holding f _G+C constant, favors the b state, thereby increasing f_b and the amplitudes of different spectrum at 287 and 295 nm with increasing T. As noted by SG, this correlation also is not consistent across all DNAs examined, and the division by f _G+C produces absurdly large quotients for DNAs with very small f _G+C, diverging (to ∞) for any DNA consisting entirely of A⋅T bp with one or more TpT steps. Nonetheless, the correlations reported by SG suggest that clustering of bases of like kinds (i.e. either pyrimidines or purines) is of comparable or greater importance to, f _A+T = 1 − f _G+C, in shifting the a⇔b equilibrium toward the b state. I recently explored other sequence ‘cluster’ patterns that might provide useful correlations with the measured α and f_b values of several DNAs using sequence data and tools in GenBank.

Among the different properties examined, the one most strongly correlated with α and f_b was the ratio, n _alt/n _tot, of the number of bps in ‘alternating’ runs, (yr)_{n /2} and (ry)_{n /2,} with n ≥6 bp, to the total number of bp, n ^tot, in both ‘alternating’ runs and ‘like’ runs, y_n and r_n, with n ≥6 bp. Formally, n ^alt = Σ( n ≥6) n × (ν((yr)_{n /2)} + ν((ry)_{n /2})), and n ^like = Σ( n ≥6) n × (ν(y_n) + ν(r_n)), where Σ( n ≥6) denotes a sum over all n ≥ 6, and ν is the number of times that the indicated run occurs in the sequence of the 5′–3′ strand of the DNA examined, and n ^tot = n ^alt + n ^like. These ratios are presented in Table 5 along with the values of (n ^tot/N)_n ≥ 6, α, and f_b. The columns are positioned in the order of monotonically increasing values of α. The end-to-end anti-correlation between the mostly declining values of (n ^alt/n ^tot) and the steadily increasing values of α is reasonably good, with two notable exceptions, pUC8 and (GC)_N⋅(GC)_N, to be discussed below. The evident correlation suggests that, for n ≥6, ‘alternating’ runs favor the a state of lower α, and like runs favor the b state of higher α. It was found that if n ^alt and n ^like include runs with all possible values of n ≥2 (instead of n ≥6), then almost no difference between ϕ29 and pBR322 DNAs is found in the n ^alt/n ^tot ratio (data not shown). The ratio (n ^tot/N)_n ≥ 6, must also be important, as it determines the mean spacing between runs. These findings suggest that sequence preferences might arise primarily from fewer, longer, and more isolated, but not too isolated, runs that significantly favor (stabilize) either the a or b states, and influence the behavior of their surrounding DNA via the large cooperativity of the two-state transition.

Table 5. Comparison of sequence properties of different DNAs with their measured α and f_b

The n ^alt are the numbers of bps in ‘alternating’, (yr)_n plus (ry)_n, runs and the n ^tot are the total numbers of bps in both ‘alternating’ and ‘like’, y_n plus r_n, runs, all with n ≥ 6 bp. The α and associated f_b values come from FPA measurements in 0.1 and 0.01 M NaCl at 293 ± 1 K (Fujimoto and Schurr, Reference Fujimoto and Schurr1990; Wu et al., Reference Wu, Song, Clendenning, Fujimoto, Benight and Schurr1988; Song et al., Reference Song, Fujimoto, Wu, Thomas, Shibata and Schurr1990; Wu et al., 1991; Delrow et al., 1998) and were corrected to apply for a fast-relaxing persistence length of 170 nm in 0.1 M NaCl (Okonogi et al., 1999; Naimushin et al. 2000). The DNAs and their associated columns are ordered in the direction of increasing α. The * on the α and f_b for M13mp7 DNA indicate that the recently linearized DNA might not have equilibrated sufficiently after relaxation of superhelical stress (Song et al. Reference Song, Fujimoto, Wu, Thomas, Shibata and Schurr1990; Brewood et al., Reference Brewood, Delrow and Schurr2010).

The n ^alt/n ^tot ratio of pUC8 should place it to the right of both control DNA and pBR322 in Table 5. Its appearance on their left side is likely due to its significantly greater (n ^tot/N)_n ≥ 6 ratio, which means that a greater fraction of all bps in the DNA are involved in n ≥6 runs of both kinds than is the case for pBR322 or control DNA. The average size of n ≥6 runs lies in the range of 6.3–7.4 bp for practically all of the DNAs examined (data not shown). The average number of bp between runs with n ≥ 6 is found to be <(7.4)((1/(n ^tot/N)_n ≥ 6) − 1) = 30, 57, and 39 bp for, respectively, pUC8, control, and pBR322 DNAs, all well within estimates of the domain sizes for the two-state model.

The (GC)_N⋅(GC)_N DNA is a special case. The correlation ‘rule’ predicts that it should exhibit a smaller α and f_b than any other DNA. Although its α and f_b exceed the corresponding values of ϕ29 DNA, which is largely in the a state, they are still smaller than those of all the other DNAs. The reported persistence length of (GC)_N⋅(GC)_N, P = 83 nm, exceeds the usual range, P = 50 ± 4 nm (Thomas and Bloomfield, Reference Thomas and Bloomfield1983), as would be expected for a DNA largely in the a state, but it also exceeds the maximum value, P = 69 nm, attainable by the two-state model. Perhaps the elastic constants of the a and b states increase somewhat with a sufficient increase in f _G+C. In any case, the suggested ‘rules’ may not apply to synthetic sequences that deviate greatly from those of most natural DNAs.

In summary, overall sequence preferences for the a and b states might be determined primarily by n ^alt/n ^tot and (n ^tot/N)_n ≥ 6, which are determined simply by the relative numbers of bps in the ‘alternating’ and ‘like’ runs with n ≥ 6, and also by the fraction of all bps in the sequence that is contained in such runs. A modest contribution of decreasing ionic strength to favor the longer b state is also suggested by, and consistent with, the available data. Secondarily, increasing values of f _G+C may modestly favor f_a, and also raise somewhat both torsion and bending elastic constants of both a and b states. These preference guidelines are adopted as the tentative working hypothesis in subsequent sections.

The present inference that ‘alternating’ runs with n ≥6 favor the a state of lower α and higher κ_β, whereas ‘like’ runs favor the b state of higher α and lower κ_β, apparently conflicts with the findings of Okonogi et al. (Reference Okonogi, Alley, Reese, Hopkins and Robinson2002), namely that ‘like’ runs with n ≥7 bp exhibit a higher κ_βf, and ‘alternating’ runs exhibit a lower κ_βf. One possible explanation for this difference is that the κ_β values of the present analysis are equilibrium values, which include permanent and slowly relaxing bends, as well as fast-relaxing dynamic bends, whereas the κ_βf values of Okonogi et al. pertain solely to the fast-relaxing dynamic bends at very short times in their experiments. The fast-relaxing κ_βf values are much greater than the equilibrium values, and contribute much less to the total mean-squared curvature, so they easily could be smaller in the a state and greater in the b state, as observed by Okonogi et al. while at the same time the equilibrium κ_β values could be smaller in the b state and greater in the a state, as found previously (Schurr, Reference Schurr2019a). Another possible explanation is that the DNAs bearing both kinds of run were assumed to have the same contour length in the calculation of the rotational relaxation times for end-over-end tumbling, which were used in the fitting to determine the elastic constants associated with the fast-relaxing bends. If the epr data of Okonogi et al. were reanalyzed by taking into account the different contour lengths of the DNAs with each kind of run, which imposes its preferred a or b state on the flanking DNA, that could conceivably reverse the conclusions regarding which kinds of runs are associated with the larger and smaller fast-relaxing bending elastic constants.

A possible explanation for the low and T-independent values of α for ϕ29 DNA

The torsion elastic constant of ϕ29 DNA was measured to be α = (4.55 ± 0.25) × 10⁻¹⁹ J in 0.1 M NaCl at 293 K, and this value remained essentially constant within experimental errors from 283 to 333 K (Thomas and Schurr, Reference Thomas and Schurr1983). Inserting the measured value of α into Eq. (8) yields, f_b = 0.0378 and f_a = 0.962. Because the measured value is so close to α_a, small errors in either α or α_a can cause very large relative errors in f_b, but not in f_a. In any event, this DNA is almost entirely in the a state at 293 K, so little further variation of its torsion or bending elastic constant with increasing T is expected, as observed. In accord with this prediction, DLS experiments on ϕ29 DNA over a wide range of scattering vectors (together with the constant α values) demonstrated that both the equilibrium bending rigidity and dynamic bending rigidity varied at most only slightly upward, from 273 to 343 K (Wilcoxon and Schurr, Reference Wilcoxon and Schurr1983).

Anti-correlation of torsion and bending elastic constants limits variations in E _T

Equation (11) yields, α = 6.63 × 10⁻¹⁹ J, for an unstrained plasmid DNA in ~40 mM univalent cations plus 5.5 mM Mg²⁺ at 310 K, which is similar to buffers used to measure topoisomer distributions and E _T values of various DNAs. Suppose that DNA with a different sequence in the same buffer could exhibit a value as low as 0.82 × 6.63 × 10⁻¹⁹ = 5.4 × 10⁻¹⁹ J. Then the corresponding values of E _T for a 4363 bp DNA reckoned for these two values of α are E _T = 1140 for α = 6.63 × 10⁻¹⁹ J and E _T = 1111 for α = 5.4 × 10⁻¹⁹ J. Thus, a 0.82-fold decrease in α produces only a 0.975-fold decrease in E _T. Evidently the decrease in α is nearly compensated by the corresponding increase in κ_β and decrease in h that accompanies the shift in the a⇔b equilibrium due to the change in sequence. Thus, DNAs with a wide range of sequences and rather different values of α and κ_β are expected to yield surprisingly similar values of E _T.

j-factors obtained from cyclization rates

Measurements of the so-called j-factor have often been used to assess helix repeats (bp/turn), persistence lengths, non-relaxing (permanent) bends, flexibilities of certain subsequences, and other properties of DNA (Shore et al., Reference Shore, Langowski and Baldwin1981; Shore and Baldwin, Reference Shore and Baldwin1983; Taylor and Hagerman, Reference Taylor and Hagerman1990; Vologodskaia and Vologodskii, Reference Vologodskaia and Vologodskii2002; Du et al., Reference Du, Smith, Shiffeldrim, Vologodskaia and Vologodskii2005; Geggier and Vologodskii, Reference Geggier and Vologodskii2010; Geggier et al., Reference Geggier, Kotlyar and Vologodskii2011). The j-factor is the concentration of the sticky end at the tail of a linear monomer with the correct position and orientation to attach to the complementary sticky end at the head of that same monomer. The j-factor is typically assumed to be the equilibrium constant for monomers to form non-covalently attached (NCA) circular monomers divided by that for monomers to form NCA linear dimers, after accounting for the three distinct kinds of dimers and their respective topological symmetry numbers (Gilson and Irikura, Reference Gilson and Irikura2010). Because of significant strain in small NCA monomer circles, which is absent in the linear species, this assumption cannot be strictly valid. Nevertheless, the prevailing equilibria in solution are sampled by a DNA ligase that seals one of the two nicks in each NCA species. The equilibrium ratio of NCA circular monomers to NCA linear dimers in solution is determined from the ratio of their ligated products extrapolated to zero reaction time. Unfortunately, j-factors obtained by this cyclization rate method rest upon one questionable assumption, and theoretical j-factors rest upon another, as described below. Hence, a comparison of theoretical and experimental j-factors does not reliably test the cooperative two-state model.

The actual j-factor, j _act, is related to the experimental value by, j _act = R × j _exp, wherein, j _exp = 2C(0) lim(t → 0) P ₁(t)/P ₂(t), is the experimental value, C(0) is the concentration of ligatable monomers at time t = 0, and P ₁(t)/P ₂(t) is the accumulated product ratio (Shore and Baldwin, Reference Shore and Baldwin1983) (SB). The factor R consists of ratios of kinetic and Michaelis constants pertaining to complexes of the ligase with NCA monomer circles on the one hand and NCA linear dimers on the other hand (Schurr, Reference Schurr2019a; Supplementary Information S5). R is assumed to be 1.0 for these doubly nicked species, but is that really true?

Measured initial rates of ligation of 247, 250, and 1361 bp circular DNAs containing "only one or a few" nicks to produce covalently closed circles were fairly similar, although the apparent first-order rate constant for the more torsionally strained 247 bp DNA was about 0.87-fold smaller than that of the other two DNAs (Reference Shore and BaldwinSB). Such similarity was taken as evidence that bending and/or torsional strains did not significantly influence the ligation rate of singly nicked circular DNAs. The steeper and more varied slopes of the ln[C(t)/C(0)] versus t plots at long times were not addressed. An analysis of the kinetics of ligation of a hypothetical mixture of DNAs containing both singly and multiply nicked circular species revealed that, at the earliest times, only that fraction of all DNAs that were initially singly nicked (i.e. not multiply nicked) contributed to the initial slope, resulting in an apparent rate constant that was the true overall first-order rate constant multiplied by the fraction (<1.0) of initially singly nicked DNAs (unpublished calculations). At sufficiently long times, all of the surviving ligatable DNAs were singly nicked, and their slopes reflected the true first-order rate constants. At long times, the 250 bp and 1361 bp samples exhibited significantly higher rate constants than the more torsionally strained 247 bp DNA. The data also hint that the 1361 bp DNA is ligated at a higher rate than the 250 bp DNA at long times, but noise in the data renders such a conclusion equivocal. The data of SB apparently support the notion that increasing strains produce significant reductions in the true first-order ligation rate constants of singly nicked circular monomers. Even larger reductions would be expected for the doubly nicked NCA species.

A detailed analysis of the likely effects of tension on the reaction kinetics suggests that for ~ 200 bp DNAs, R ≥1.0, and that upon increasing T, R should decline toward 1.0, and j _act should approach j _exp (Schurr, Reference Schurr2109a, Supporting Information S5).

The effects of kinks, which are likely promoted by tension in small NCA circular monomers, which are doubly nicked, would act to increase the j _act above the theoretical value, j _thc, reckoned for a DNA with the correct elastic constants but a circular (as opposed to kinked) ground state. Moreover, with decreasing T, the probability of kink formation should decline and j _act should approach j _thc.

The inequalities, j _act ≥ j _exp and j _act ≥ j _thc, of the two preceding paragraphs provide no indication of the relative magnitudes of j _thc and j _exp. However, the first inequality is expected to approach its equality at higher T, where R approaches 1.0, whereas the latter inequality is expected to approach its equality at lower T, where the kink probability is lower. Hence, it is possible that j _thc could significantly exceed j _exp at a lower T and that j _exp could exceed j _thc at higher T.

Comparison of the two-state model with measured j-factors

Using an appropriate theory for DNAs with circular ground states (Shimada and Yamakawa, Reference Shimada and Yamakawa1984), values of j _thc for the present two-state model at various temperatures were reckoned for DNAs with lengths chosen to have precisely 19.0 turns at each prevailing T (Schurr, Reference Schurr2019a). These j _thc values were compared with j _exp values that were calculated for molecules of the same lengths with the temperature-dependent persistence lengths reported by Geggier et al. by using the same torsion elastic constant, α = 9.12 × 10⁻¹⁹ J, and rise per bp, 0.34 nm, that they used to extract such persistence lengths from their measured j _exp values (cf. Table 6 of Schurr, Reference Schurr2019a). As anticipated above, j _thc exceeded j _exp (by 2.8-fold) at 278 K, whereas j _exp exceeded j _thc (by 1.9-fold) at 315 K, and they matched between 303 and 310 K. Subsequently, I recalculated j _thc values for the same DNAs and temperatures using the same two-state model, but by evaluating the exact partition function for small circles and taking full account of the lower eigenvalue. In a comparison against the same experimental values, the new j _thc exceeds j _exp (by 1.8-fold) at 278 K, whereas j _exp exceeds the new j _thc (by 2.2-fold) at 315 K, and they match just below 295 K. Are the discrepancies at 278 and 315 K due to a failure of the two-state model, or instead to a failure of assumptions employed in the measurement and calculation of j-factors for the doubly nicked circular DNAs? Because the two-state model gives such a good account of the torsion elastic constant and supercoiling free energies of other small circular DNAs that are free of nicks at 293 and 310 K, it seems likely that such disagreements are attributable in large part to unequal kinetics constants of the doubly nicked monomer circles and linear dimers and/or to a kinked or otherwise non-circular ground state of doubly nicked monomer ‘circles.’

j-factors from single-molecule fluorescence resonance energy transfer (FRET)

Estimation of j-factors from single-molecule FRET measurements is also problematic (Schurr, Reference Schurr2019a, Supporting Information S5). The difficulty with the circular ground state in the theory for NCA monomer circles remains. Typically, in FRET, the ratio of equilibrium constants to form NCA monomers from linear monomers, to that for forming short NCA linear dimers from short linear monomers is reckoned from the ratio of their annealing rates. This rate method is confronted by new problems rendering quantitative interpretation uncertain (Mulligan et al., Reference Mulligan, Chen, Phillips and Spakowicz2015; Jeong et al., Reference Jeong, Le and Kim2016). Also the nearby surface may differently affect the equilibrium constants for formation of annealed monomers and annealed linear dimers of short DNAs.

In summary, both j-factor measurements and persistence lengths obtained therefrom are quantitatively suspect in regard to absolute values but might still provide useful information regarding relative values for different sequences and conditions.

Effects of a small change in a much larger sequence

The discovery that a small change in a much larger DNA sequence could produce disproportionately large effects on its properties came from an experiment originally intended to examine the flexibility associated with B–Z junctions (Kim et al., Reference Kim, Fujimoto, Furlong, Sundstrom, Humbert, Teller and Schurr1993). Twenty-five consecutive bps near the middle of the ~1100 bp HincII–HincII restriction fragment of pBR322 were replaced by 16 bp of alternating GC, or (GC)₈. Less than 2.3% of the original sequence of the DNA was altered. Subsequent sequencing of the resulting 1097 bp insert DNA revealed that the final C of the (GC)₈ had been dropped during the construction/amplification process, leaving a 16 bp run of alternating CG, the first C of which was part of the original flanking sequence. The control DNA contained the unaltered sequence, but was 8 bp shorter, due to clipping with HinfI at its 3′-end to facilitate purification. Both insert and control DNA samples had no detectable levels of protein and/or polyamine contaminants in fluorescamine tests, exhibited no detectable aggregation in DLS or sedimentation velocity measurements, and after completion of all experiments exhibited no detectable levels of single-strand breaks in denaturing gels. Both insert and control DNAs were normal in regard to (i) UV spectra, (ii) electrophoretic mobilities compared with ladders of two different lengths, (iii) sedimentation coefficients, and (iv) translational diffusion coefficients. The apparent DLS diffusion coefficients of the two samples were very similar at all scattering vectors (K) from zero up to the plateau region at K ² = 20 × 10⁻¹⁰ cm⁻². However, in other regards, the insert and control samples were strikingly different, even in 0.1 M NaCl at ambient T (293–294 K). The CD band of the insert DNA near 273 nm exceeded that of the control DNA by ~1.3-fold, suggesting a large increase in f_a caused by the insert. The α-values of the 1089 bp control DNA were highly non-uniform over the different (20–120 ns) time scales measured by FPA, indicative of a large domain of substantially higher α than the effective value (α ~ 6.0 × 10⁻¹⁹ J) prevailing at long times, t ≥ 80 ns. In contrast, the torsion elastic constant of the insert DNA (~4.55 × 10⁻¹⁹ J) was much more uniform, and for t ≥ 80 ns was smaller than that of the control DNA by a factor of ~1.32–1.35, which is further evidence of a large increase in f_a. The insert DNA was also more susceptible than the control DNA to cleavage by the random single-strand cutter, S1 nuclease, over almost its entire length.

The two-state cooperative model was applied to better understand these disproportionately large effects of the sequence change. Additional assumptions regarding the preferences of subsequences of these DNAs for the a and b states were required. The 16 bp alternating CG insert had a large fraction of G⋅C bp, f _G+C = 1.0, and comprised a rather long ‘alternating’ (yr)₈ run, both of which should substantially favor the a state. It was assumed to exhibit always the a state.

The original 25 bp sequence in the control DNA is characterized by f _G+C = 0.375, which should modestly favor the b state, and one y₃ run partially overlapped by one (yr)₂ run, but no longer runs of either kind. However, an r₅ run terminates just 4 bp upstream from the 5′-end of that 25 bp sequence, and a y₆ run terminates several bps upstream from that. Thus, the 25 bp sequence is proximal to a region with a rather strong preference for the b state and is assumed to always adopt the b state. The 16 bp alternating (CG)₈ sequence that replaces the 25 bp sequence is assumed to dominate the nearby ‘like’ runs, so as to effectively replace the initial domain of b in the control DNA with a 16 bp domain of a in the insert DNA.

For simplicity, the 25 bp control and 16 bp insert sequences were assumed to be located at the precise centers of DNAs of the same length as their actual lengths, so only half-DNAs with appropriate boundary conditions at one end (near the middle of the full molecule) needed to be considered. The cooperativity parameter, J, in the M ₁₂ position of the transfer matrix, M, was replaced by, J_ab, and that in the M ₂₁ position by, J_ba. Both J_ab and J_ba were eventually assigned the same value as J in Table 1, but this artifice enabled determination of the average numbers of ab and ba junctions separately, from which average sizes of a and b domains were approximately inferred. The flanking DNAs were assumed to have a uniform value of B, which was chosen to provide the value, α = 5.01 × 10⁻¹⁹ J, when such DNA was not connected to the central subsequences. This choice allowed the best overall agreement with the experimental data.

The acting torsion elastic constant was simply calculated from f_b, α_b, and α_b, and is denoted by α_{u 1}. An equal number of torsion springs, each with uniform elastic constant, α_{u 1}, would exhibit the same end-to-end torsion elastic constant as the original chain of different springs. ‘Average’ distributions of the a and b states for control and insert DNAs are shown in Fig. 5.

Fig. 5. Schematic illustration of the positions and average sizes of the different domains of a (light gray) and b (dark gray) states in the control (top) and insert (bottom) DNAs that are predicted by the two-state cooperative model.

The numbers of springs, average domain sizes of the central domains, and of the second and third outer domains for both control and insert DNAs, and their α_{u 1} values are listed in Table 6. The small insert has effectively switched a 252 bp central domain of state b to state a, thereby reducing the b content of the DNA to slightly less than half its original value. This results in a reduction of α_{u 1} from 5.70 to 4.98, as indicated in Table 6. This is a surprisingly large decrease for such a small change in sequence, although it is only about half as large as the decrease observed in the FPA.

Table 6. Calculated properties of different DNA models in 0.1 M NaCl at 293 K

The Nspr are the numbers of springs between the subunits of each model. The n ₁ are the mean number of bp in the central domain, which is state b in control DNA, and is state a in insert DNA. The n ₂ and n ₃ are the numbers of subunits in the second and third domains outward from the center (counting those on both sides). The α_{u 1} and α_{u 2} values are reckoned as described previously (Schurr, Reference Schurr2019b).

A second kind of uniform torsion elastic, α_{u 2}, is the value required for a chain of uniform springs to produce the same value of the FPA as the actual model DNA at a particular time, t = 100 ns, after the polarized exciting pulse. This time was chosen because extension of the time span of the experimental data analyzed from 0 to 80 ns to 0 to 120 ns typically produces no significant change in the best-fit value of α. The FPA and α_{u 2} are sensitive not only to the fraction of the bp in the b state, but also to the distribution of stiffer b and softer a springs along the DNA. In general, α_{u 2} differs at least slightly from α_{u 1}. The computational effort required to calculate the mean FPA decay at any single time for ~1100 bp DNAs for all of the different appropriately weighted configurations of the a and b states would be extreme, requiring the determination of eigenvalues and eigenvectors for a colossal number of ~1100 × 1100 torque constant matrices. Hence, a far simpler and more approximate approach was adopted (Schurr, Reference Schurr2019b). Specifically, the length of the DNA was drastically shortened to 23 subunits separated by 22 bending springs. The length of each subunit was scaled up so as to retain the length of the original DNA, and the values of the elastic constants of the springs between bps were correspondingly scaled down, so as to provide the correct end-to-end torsional elastic constant of the entire DNA. The hydrodynamic cylinder volume per subunit, which governs the friction coefficient for azimuthal rotation, was assumed to be the same for both a and b states. These reduced models, called Redcontrol or Redinsert, exhibited f_b values and relative domain sizes comparable to those of the full-length DNAs. The torsion elastic constants of the springs at the ab and ba junctions were taken to be α_a.

The FPA at time, t = 100 ns, after the exciting pulse was reckoned for each of two closely related ‘typical’ distributions of states among the 23 subunits of both Redcontrol and Redinsert models by using all of the eigenvalues of the particular 22 × 22 torque constant matrix in each case (Schurr, Reference Schurr2019b). Then α_{u 2} for the corresponding uniform torque constant matrix that gave the same FPA at t = 100 ns was determined by successive adjustment to achieve a matching value of the FPA at t = 100 ns. These values are presented in Table 6.

For all reduced models, α_{u 2} exceeds α_{u 1}. More importantly, the α_{u 2} values of the two Redcontrol models exceed the larger Redinsert α_{u 2} value (with f_b = 0.182) by factors of 1.21 and 1.22, and exceed the smaller of the Redinsert α_{u 2} values (with f_b = 0.136) by factors of 1.26 and 1.27. Though still below the ratio of measured values for the control and insert DNAs, this implies that the two-state model is in both qualitative and approximately quantitative accord with the experimental results. If J is assumed to be half as large as the value in Table 1, the calculated ratio of α_{u 2} values for insert and control DNAs agrees very well with the experimental ratio, ~1.35 (not shown). Given the likelihood that J is also sequence dependent, this possibility cannot be ruled out.

Given the substantial differences in other properties, the similarity of the apparent diffusion coefficients of the insert and control DNAs is puzzling. From the fractions of bps in the b state, the acting bending rigidities and persistence lengths of the two DNAs were calculated, and the ratio of their diffusion coefficients at K ² = 0 was computed, as described previously (Yamakawa and Fujii, Reference Yamakawa and Fujii1973; Schurr, Reference Schurr2019a). The predicted ratio, D ₁₀₉₇/D ₁₀₈₉ = 0.981, matches the experimental results well within their experimental errors. The close agreement between the diffusion coefficients of the two DNAs is attributable in part to the free-draining character of such short DNAs and in part to the fact that the modestly greater persistence length and slightly larger number of bp of the insert DNA are significantly offset by its smaller mean rise per bp.

In summary, the two-state cooperative model with certain assumptions about sequence preferences for the a and b states gives a good qualitative account of the switching of a ~252 bp domain from state b in the control DNA to a in the insert DNA, and the reduced models provide an approximately quantitative estimate of the ratio of α_{u 2} values for the control and insert DNAs.

Discrepant measurements of intrinsic curvature

Persistence length(s) of synthetic straight and natural DNAs were measured by cryoEM of samples at 101 K (Bednar et al., Reference Bednar, Furrer, Katritch, Stasiak, Dubochet and Stasiak1995) and later by j-factor measurements at 295 K (Vologodskaia and Vologodskii, Reference Vologodskaia and Vologodskii2002). The synthetic straight DNAs comprised tandem repeats of a 5 bp sequence, so any largely planar bend in the first 5 bp of a turn was compensated by that in the second 5 bp. The vast majority of natural DNAs examined in the cryoEM experiment and those studied in the j-factor experiments contained no strong bending sequences, so those experiments concerned the proposed small sequence-dependent permanent bends at each of the 16 different bp steps (Bolshoy et al., Reference Bolshoy, McNamara, Harrington and Trifonov1991). The cryoEM experiments found much larger persistence lengths for straight DNAs than for natural DNAs, indicating that permanent bends (intrinsic curvature) made a very substantial contribution to the mean-squared bending angle(s) of natural DNAs, which are proportional to their inverse persistence lengths. In contrast, the j-factor results seemed to indicate only a very small contribution of permanent bends to the mean-squared bending angle(s) of natural DNAs. This large discrepancy was resolved by considering the different conditions prevailing in the two experiments and imposing additional assumptions on the two-state model (Schurr, Reference Schurr2019b).

The value, α = 5.6 × 10⁻¹⁹ J, measured in 4 mM univalent cations (Naimushin et al., Reference Naimushin, Fujimoto and Schurr2000) was assumed to apply in the cryoEM buffer, which contained ~5 mM univalent cations. From that value, the f_b, κ_β, and the total persistence length, P _tot, of the natural DNA were reckoned using the two-state model parameters. The persistence length associated with the fast-relaxing bends was taken to be the measured value, P_f = 200 nm, in ~4 mM univalent cations (Naimushin et al., Reference Naimushin, Fujimoto and Schurr2000), and was assumed to apply to both a and b states. It was assumed that the a state exhibited non-relaxing (permanent) bends, but no slowly relaxing bends, and that the b state exhibited slowly relaxing bends, but no non-relaxing bends. These assumptions enabled a determination of the fast-relaxing, slowly relaxing, and non-relaxing components of the persistence lengths (P_a and P_b) of the individual a and b states of the natural DNA, from which its P _tot could be reckoned, although in that case it was already determined from the assumed value of α. The two-state model, augmented by these assumptions, yielded a persistence length associated with non-relaxing bends, P_n = 134⋅6, close to the generic value, P_n = 137 nm, estimated from the proposed bp step values (Schellman and Harvey, Reference Schellman and Harvey1995). The straight DNA was assumed to exhibit the same values of α, f_b, and P_b as the natural DNA, but its P_a = 200 nm. These assumptions enabled appropriate estimates of the quantities in the initial solution at 293 K. The drop in T upon immersion in the cryogenic fluid was assumed to be so rapid that the a⇔b equilibrium did not shift, and the slowly relaxing bends did not relax. Thus, persistence lengths associated with the slowly relaxing bends, as well as the non-relaxing bends, did not change. However, some motion is allowed down to 223 K, below which water rapidly freezes or vitrifies, so equilibration of the fast-relaxing bending motions down to 223 K was assumed to occur. The persistence length, P_f = 200 nm, associated with the fast-relaxing bends at 293 was scaled by the ratio 293/223 to apply at 223 K. A comparison of the theoretical and experimental values of P _tot for the synthetic straight and natural DNAs in the cryoEM experiment is presented in Table 7.

Table 7. Comparison of theoretical (P ^th) and experimental (P ^exp) persistence lengths for straight and natural DNAs under the different conditions of cryoEM and j-factor measurements

It was assumed that non-relaxing bends occur only in the a state and slowly relaxing bends occur only in the b state.

Although the theoretical cryoEM values for both natural and straight DNAs significantly exceed their corresponding experimental values, the theoretical ratio (101.6/57.6 = 1.75) is close to the experimental ratio (82/46 = 1.78). Some of the discrepancy between theoretical and measured values for each DNA might be attributed to an unexpected decline of the fast-relaxing bending elastic constants of the a and b states with decreasing temperature below 278 K in the experiments, since 223 K is far below the range for which such ‘constant’ values were determined.

The j-factor measurements were made in ~40 mM univalent cations plus 10 mM Mg²⁺ at 295 K. The linear forms of both natural and synthetic straight DNAs in this buffer were assumed to exhibit, α ₀ = 7.55 × 10⁻¹⁹ J, from Eq. (11). From this value, the f_{b 0}, κ_{β 0}, and P _tot of the linear natural DNA were reckoned. The a state of the linear straight DNA was assumed to contain neither non-relaxing nor slowly relaxing bends. In this buffer, its fast-relaxing bends were assumed to yield, P_a = 150 nm, but its P_b remained the same as for the natural DNA in the same buffer. The P _tot for the linear straight DNA could be reckoned directly from its f_b, P_a, and P_b, the latter two of which were assumed to remain constant.

Theoretical j-factors were reckoned beginning with, α ₀ = 7.55 × 10⁻¹⁹ J, following the general approach described previously (Schurr, Reference Schurr2019b), but now using the full circle boundary conditions and taking account of the lower eigenvalue, as discussed in sections ‘Development of the two-state model’ and ‘Tests of the cooperative two-state model vs. experiment'. Persistence lengths were computed from the theoretical j-factors using the same protocol and assumptions as Vologodskaia and Vologodskii (Reference Vologodskaia and Vologodskii2002). These theoretical persistence lengths for the straight and natural DNAs are compared with the experimental values in Table 7. Agreement of the theoretical and experimental persistence lengths is remarkably good, although the significance of that is uncertain, given the previously noted problems with both experimental and theoretical j-factors.

The good agreement between the ratio, P ^th(strt)/P ^th(nat) = 1.05, where ‘strt’ denotes ‘straight’ and ‘nat’ denotes ‘natural’, and P ^exp(strt)/P ^exp(nat) = 1.03 strongly suggests that the theory has identified the reason for the near absence of permanent bends in the j-factor experiments. That behavior is due to the enhanced fraction, f_b = 0.868, of the bp in the b state, which arises from the different buffer and the tension in the small circles, and also to the assumed absence of permanent bends in the b state. This latter assumption was predicated on the notion that the bps are less restricted in the a state, as evidenced by its higher entropy. That might allow them to deviate from the configuration preferred by the backbone toward that preferred by the bps, leading to permanent bends, whereas the same bps in the b state might have to conform more strictly to the configuration(s) dictated by the backbone, which likely suppress permanent bends. The slow relaxation of most of the mean-squared bending of the b state likely arises from a series of potential barriers along the bending coordinate(s) that are equivalent to a high internal friction.

The assumed absence of non-relaxing bends in the b state is consistent with the lack of any indication of the ‘pulling out’ of permanent bends at forces >2.0 pN in z/z _max versus F experiments.