# Statistical Mechanics of Torque Induced Denaturation of DNA

## Abstract

A unifying theory of the denaturation transition of DNA, driven by temperature $`T`$ or induced by an external mechanical torque $`\mathrm{\Gamma }`$ is presented. Our model couples the hydrogen-bond opening and the untwisting of the helicoidal molecular structure. We show that denaturation corresponds to a first-order phase transition from B-DNA to d-DNA phases and that the coexistence region is naturally parametrized by the degree of supercoiling $`\sigma `$. The denaturation free energy, the temperature dependence of the twist angle, the phase diagram in the $`T,\mathrm{\Gamma }`$ plane and isotherms in the $`\sigma ,\mathrm{\Gamma }`$ plane are calculated and show a good agreement with experimental data.

PACS Numbers : 87.14.Gg, 05.20.-y, 64.10.+h

Denaturation of the DNA, due to its essential relevance to transcription processes has been the object of intensive works in the last decades. Experiments on dilute DNA solutions have provided evidence for the existence of a thermally driven melting transition corresponding to the sudden opening of base pairs at a critical temperature $`T_m`$ . Later, following the work of Smith et al., micromanipulation techniques have been developed to study single-molecule behaviour under stress conditions and how structural transitions of DNA can be mechanically induced. While most single-molecule experiments have focused on stretching properties so far, the response of a DNA molecule to an external torsional stress has been studied very recently , sheding some new light on denaturation . From a biological point of view, torsional stress is indeed not unusual in the living cell and may strongly influence DNA functioning.

For a straight DNA molecule with fixed ends, the degree of supercoiling $`\sigma =(TwTw_0)/Tw_0`$ measures the twist $`Tw`$ (i.e. the number of times the two strands of the DNA double-helix are intertwined) with respect to its counterpart $`Tw_0`$ for an unconstrained linear molecule. In Strick et al. experiment , a $`\lambda `$ DNA molecule, in 10 mM PB, is attached at one end to a surface and pulled and rotated by a magnetic bead at the other end. At stretching forces of $`0.5`$ pN, sufficient to eliminate plectonems by keeping the molecule straight, a torque induced transition to a partially denaturated DNA is observed. Beyond a critical supercoiling $`\sigma _c0.015`$ and an associated critical torque $`\mathrm{\Gamma }_c0.05eV/\text{rad}`$ the twisted molecule separates into a pure B-DNA phase with $`\sigma =\sigma _c`$ and denaturated regions with $`\sigma =1`$. Extra turns applied to the molecule increase the relative fraction of d-DNA with respect to B-DNA.

In this letter, we provide a unifying understanding of both thermally and mechanically induced denaturation transitions. We show that denaturation can be described in the framework of first-order phase transitions with control parameters being the temperature and the external torque. This is in close analogy to the liquid-gas transition, where control parameters are the temperature and the pressure. Our theory gives a natural explanation to the BDNA-dDNA phases coexistence observed in single molecule experiments . We give quantitative estimates for the denaturation free-energy $`\mathrm{\Delta }G`$, the temperature dependence of the average twist angle $`\mathrm{\Delta }\theta /\mathrm{\Delta }T`$, the critical supercoiling $`\sigma _c`$ and torque $`\mathrm{\Gamma }_c`$ at room temperature in good agreement with the experimental data. Furthermore the dependence of the critical torque as a function of the temperature is predicted.

Our model reproduces the Watson-Crick double helix (B-DNA) as schematized fig 1. For each base pair ($`n=1,\mathrm{},N`$), we consider a polar coordinate system in the plane perpendicular to the helical axis and introduce the radius $`r_n`$ and the angle $`\phi _n`$ of the base pair . The sugar phosphate backbone is made of rigid rods, the distance between adjacent bases on the same strand being fixed to $`L=6.9\AA `$. The distance $`h_n`$ between base planes $`n1`$ and $`n`$ is expressed in terms of the radii $`r_{n1},r_n`$ and the twist angle $`\theta _n=\phi _n\phi _{n1}`$ as

$$h_n(r_n,r_{n1},\theta _n)=\sqrt{L^2r_n^2r_{n1}^2+2r_nr_{n1}\mathrm{cos}\theta _n}.$$
(1)   
The potential energy associated to a configuration of the degrees of freedom $`(r_n,\phi _n)`$ is the sum of the following nearest neighbor interactions.

First, hydrogen bonds inside a given pair $`n`$ are taken into account through the short-range Morse potential  $`V_m(r_n)=D\left(e^{a(r_nR)}1\right)^2`$ with $`R=10\AA `$ . Fixing $`a=6.3\AA ^1`$ , the width of the well amounts to $`3a^10.5\AA `$ in agreement with the order of magnitude of the relative motion of the hydrogen bonded bases . A base pair with diameter $`r>r_d=R+6/a`$ may be considered as open. The potential depth $`D`$, typically of the order of $`0.1eV`$  depends on the base pair type (Adenine-Thymine (AT) or Guanine-Citosine (GC)) as well as on the ionic strength.

Secondly, the shear force that opposes sliding motion of one base over another in the B-DNA conformation is accounted for by the stacking potential  $`V_s(r_n,r_{n1})=Ee^{b(r_n+r_{n1}2R)}(r_nr_{n1})^2`$. Due to the decrease of molecular packing with base pair opening, the shear prefactor is exponentially attenuated and becomes negligible beyond a distance $`5b^1=10\AA `$, which coincides with the diameter of a base pair .

Thirdly, an elastic energy $`V_b(r_n,r_{n1},\theta _n)=K[h_nH]^2`$ is introduced to describe the vibrations of the molecule in the B phase. The helicoidal structure arises from $`H<L`$: in the rest configuration $`r_n=R`$ at $`T=0`$K, $`V_b`$ is minimum and zero for the twist angle $`\theta _n=2\pi /10`$. Choosing $`H=3\AA `$, we recover at room temperature $`T=298`$K the thermal averages $`h_n3.4\AA `$ and $`\theta _n2\pi /10.4`$ . The above definition of $`V_b`$ holds as long as the argument of the square root in (1) is positive, that is if $`r_n,r_{n1},\theta _n`$ are compatible with rigid rods having length $`L`$. By imposing $`V_b=\mathrm{}`$ for negative arguments, unphysical values of $`r_n,r_{n1},\theta _n`$ are excluded. As the behaviour of a single strand ($`r>r_d`$) is uniquely governed by this rigid rod condition, the model does not only describe vibrations of helicoidal B-DNA but is also appropriate for the description of the denaturated phase.

As will be discussed later, the elastic constant $`K=0.014eV/\AA ^2`$ is determined to give back the torsional modulus $`C`$ of B-DNA estimated to $`C=860\pm 100\AA `$  at $`T=298`$K. The parameters of the Morse potential $`D`$ and of the stacking interaction $`E`$ we have set to fit the melting temperature $`T_m=350`$K of the homogeneous Poly(dGdT)-Poly(dAdC)-DNA at $`20mMNa^+`$ , see inset of fig 3. This melting temperature coincides with the expected denaturation temperature of a heterogeneous DNA with a sequence GC/AT ratio equal to unity at $`10mMNa^+`$ , as the $`\lambda `$-DNA in the experimental conditions of . Among all possible pairs of parameters $`(D,E)`$ that correctly fit $`T_m`$, we have selected the pair $`(D=0.16eV,E=4eV/\AA ^2)`$ giving the largest prediction for $`\mathrm{\Delta }G`$, see inset of fig 2, that is in closest agreement with thermodynamical estimates of the denaturation free-energy.

When the molecule is fixed at one end and subject to a torque $`\mathrm{\Gamma }`$ on the other extremity, an external potential $`V_\mathrm{\Gamma }(\theta _n)=\mathrm{\Gamma }\theta _n`$ has to be included. A torque $`\mathrm{\Gamma }>0`$ overtwists the molecule, while $`\mathrm{\Gamma }<0`$ undertwists it.

The configurational partition function at inverse temperature $`\beta `$ can be calculated using the transfer integral method :

$$Z_\mathrm{\Gamma }=_{\mathrm{}}^{\mathrm{}}𝑑\phi _NR,\phi _N|T^N|R,0$$
(2)   
As in the experimental conditions, the radii of the first and last base pairs are fixed to $`r_1=r_N=R`$. The angle of the fixed extremity of the molecule is set to $`\phi _1=0`$ with no restriction whereas the last one $`\phi _N`$ is not constrained. The transfer operator entries read $`<r,\phi |T|r^{},\phi ^{}>X(r,r^{})\mathrm{exp}\{\beta (V_b(r,r^{},\theta )+V_\mathrm{\Gamma }(\theta ))\}\chi (\theta )`$ with $`X(r,r^{})=\sqrt{rr^{}}\mathrm{exp}\{\beta (V_m(r)/2+V_m(r^{})/2+V_s(r,r^{}))\}`$. The $`\sqrt{rr^{}}`$ factor in $`X`$ comes from the integration of the kinetic term; $`\chi (\theta )=1`$ if $`0\theta =\phi \phi ^{}\pi `$ and 0 otherwise to prevent any clockwise twist of the chain. At fixed $`r,r^{}`$, the angular part of the transfer matrix $`T`$ is translationally invariant in the angle variables $`\phi `$, $`\phi ^{}`$ and can be diagonalized through a Fourier transform. Thus, for each Fourier mode $`k`$ we are left with an effective transfer matrix on the radius variables $`T_k(r,r^{})=X(r,r^{})Y_k(r,r^{})`$ with

$$Y_k(r,r^{})=_0^\pi 𝑑\theta e^{\beta (V_b(r,r^{},\theta )+V_\mathrm{\Gamma }(\theta ))}e^{ik\theta }.$$
(3)   
The only mode contributing to $`Z_\mathrm{\Gamma }`$ is $`k=0`$ once $`\phi _N`$ has been integrated out in (2). The eigenvalues and eigenvectors of $`T_0`$ will be denoted by $`\lambda _q^{(\mathrm{\Gamma })}`$ and $`\psi _q^{(\mathrm{\Gamma })}(r)`$ respectively with $`\lambda _0^{(\mathrm{\Gamma })}\lambda _1^{(\mathrm{\Gamma })}\mathrm{}`$. In the $`N\mathrm{}`$ limit, the free-energy density $`f^{(\mathrm{\Gamma })}`$ does not depend on the boundary conditions on $`r_1`$ and $`r_N`$ and is simply given by $`f^{(\mathrm{\Gamma })}=k_BT\mathrm{ln}\lambda _0^{(\mathrm{\Gamma })}`$.

Note that the above result can be straightforwardly extended to the case of a molecule with a fixed twist number $`Tw=N\mathrm{}`$, e.g. for circular DNA. Indeed, the twist density $`\mathrm{}`$ and the torque $`\mathrm{\Gamma }`$ are thermodynamical conjugated variables and the free-energy at fixed twist number $`\mathrm{}`$ is the Legendre transform of $`f_\mathrm{\Gamma }`$.

We have resorted to a Gauss-Legendre quadrature for numerical integrations over the range $`r_{min}=9.7\AA <r<r_{max}`$. The Morse potential $`V_m`$ increases exponentially with decreasing $`r<R`$ and may be considered as infinite for $`r<9.7\AA `$ . The extrapolation procedure to $`r_{max}\mathrm{}`$ depends on the torque value $`\mathrm{\Gamma }`$ and will be discussed below. Using Kellog’s iterative method , the eigenvalues $`\lambda _q^{(\mathrm{\Gamma })}`$ and associated eigenvectors $`\psi _q^{(\mathrm{\Gamma })}(r)`$ have been obtained for $`q=0,1,2`$. Like a quantum mechanical wave function, $`\psi _0^{(\mathrm{\Gamma })}(r)`$ gives the probability amplitude of a base pair to be of radius $`r`$. Two quantities of interest are: the percentage of opened base pairs $`P=_{r_d}^{\mathrm{}}𝑑r|\psi _0^{(\mathrm{\Gamma })}(r)|^2`$, the averaged twist angle $`\theta =f_\mathrm{\Gamma }/\mathrm{\Gamma }`$.

Results for a freely swiveling molecule at room temperature are as follows. $`\psi _0^{(\mathrm{\Gamma }=0)}`$ is entirely confined in the Morse potential well and describes a closed molecule. Conversely the following eigenfunctions $`\psi _1^{(\mathrm{\Gamma }=0)}`$, $`\psi _2^{(\mathrm{\Gamma }=0)},\mathrm{}`$ correspond to an open molecule: they extend up to $`r_{max}`$ and vanish for $`r<r_d`$. They are indeed orthogonal to another family of excited states that are confined in the Morse potential with much lower eigenvalues. The shape of the open states are strongly reminiscent of purely diffusive eigenfunctions, $`\psi _q^{(\mathrm{\Gamma }=0)}(r)\mathrm{sin}(q\pi (rr_d)/(r_{max}r_d))`$ leading to a continuous spectrum in the limit $`r_{max}\mathrm{}`$.

This observation can be understood as follows. For $`r,r^{}>r_d`$, the transfer operator $`T_0(r,r^{})`$ is compared fig. 2 to the exact conditional probability $`\rho (r,r^{})`$ that the endpoint of a backbone rod of length $`L`$ is located at distance $`r^{}`$ from the vertical reference axis knowing that its other extremity lies at distance $`r`$ . For fixed $`r`$, $`T_0`$ and $`\rho `$ both diverge in $`r^{}=r\pm L`$ and are essentially flat in between. The flatness of $`T_0`$ derive from the expression of $`V_b`$: a rigid rod with extremities lying in $`r,r^{}`$ may always be oriented with some angle $`\theta ^{}`$ ($`0`$ at large distances) at zero energetic cost $`V_b(r,r^{},\theta ^{})=0`$. As a conclusion, our model can reproduce the purely entropic denaturated phase.

As shown fig 2, at a critical temperature $`T_m=350`$K, $`\lambda _0^{(\mathrm{\Gamma }=0)}`$ crosses the second largest eigenvalue and penetrates, as in a first-order-like transition the continuous spectrum. For $`T>T_m`$, the bound state disappears and $`\psi _1^{(\mathrm{\Gamma }=0)}`$ in fig 2 becomes the eigenmode with largest eigenvalue . The percentage of opened base pairs $`P`$ exhibits an abrupt jump from 0 to 1 at $`T_m`$, reproducing the UV absorbance vs. temperature experimental curve for Poly(dGdT)-Poly(dAdC)-DNA . The difference $`\mathrm{\Delta }G`$ between the free energy $`f_d^{(\mathrm{\Gamma }=0)}`$ of the open state ($`q=1`$) and the free energy $`f_B^{(\mathrm{\Gamma }=0)}`$ of the close state ($`q=0`$) gives the denaturation free energy at temperature $`T`$, see fig. 2. At $`T=298`$K, we obtain $`\mathrm{\Delta }G=0.022eV`$ in good agreement with the free energy of the denaturation bubble formation $`\mathrm{\Delta }G0.025eV`$ estimated in AT rich regions . The thermal fluctuations in the B-DNA phase lead to an undertwisting $`\mathrm{\Delta }\theta /\mathrm{\Delta }T\mathrm{1.4\hspace{0.17em}10}^4\text{rad/K}`$ which closely agrees with experimental measures $`\mathrm{\Delta }\theta /\mathrm{\Delta }T\mathrm{1.7\hspace{0.17em}10}^4\text{rad/K}`$ .

The presence of an overtwisting (respectively undertwisting) torque $`\mathrm{\Gamma }>0`$ (resp. $`\mathrm{\Gamma }<0`$) strongly affects $`f_B^{(\mathrm{\Gamma })}`$, leaving almost unchanged the single strand free-energy $`f_d^{(\mathrm{\Gamma })}`$. The denaturation transition takes place at $`T_m(\mathrm{\Gamma })`$ , see the phase diagram shown in the inset of fig 3. We expect a critical point at a high temperature and large positive torque such that $`\psi _1^\mathrm{\Gamma }`$ is centered on $`R`$ .

The supercoiling, induced by a torque at a given temperature smaller than $`T_m(\mathrm{\Gamma }=0)=350`$K, is the relative change of twist with respect to the value at zero torque in the B-DNA state, $`\sigma (\mathrm{\Gamma })=(\theta _\mathrm{\Gamma }\theta _{\mathrm{\Gamma }=0})/\theta _{\mathrm{\Gamma }=0}`$. In fig 3, we have plotted the isotherms in the $`\sigma ,\mathrm{\Gamma }`$ plane. Horizontals lines are critical coexistence regions between the B-DNA phase, on the left of the diagram and the denaturated phase on the right (with $`\sigma =1`$). The left steep line is found to define a linear relation between $`\mathrm{\Gamma }`$ and $`\sigma `$ : $`\mathrm{\Gamma }=K_\theta (\theta _\mathrm{\Gamma }\theta _{\mathrm{\Gamma }=0})`$. The slope $`K_\theta `$ does not vary with temperature over the range 298 K$`<T<`$350 K and is related to the torsional modulus $`C`$ of B-DNA through $`C=K_\theta h_n/(k_BT)`$. The value of $`K`$ appearing in the elastic potential $`V_b`$ and given above was tuned to ensure that $`C=860\AA `$. At room temperature, critical coexistence between B-DNA and d-DNA arises at torque $`\mathrm{\Gamma }_c=0.035eV/\text{rad}`$ and supercoiling $`\sigma _c=0.01`$. These theoretical results are in good agreement with the values $`\mathrm{\Gamma }_c=0.05eV/\text{rad}`$, $`\sigma _c=0.015`$ obtained experimentally.

We plan to combine the present model with existing elasticity theories of DNA  to understand the influence of an external stretching force on the structural transition studied in this paper. It would also be interesting to see how the above results are modified in presence of a heterogeneous sequence.

Acknowledgements : The present model is the fruit of a previous collaboration  of one of us (S.C.) with M. Barbi and M. Peyrard which we are particularly grateful to. We also thank B. Berge, D. Bensimon, C. Bouchiat, E. Bucci, A. Campa, A. Colosimo, V. Croquette, A. Giansanti for useful discussions.