Two-Phase Dynamics of DNA Supercoiling based on DNA Polymer Physics

We study the dynamics of torque propagation on supercoiled DNA based on the WLC model numerically, i.e., via an implementation of the dWLC method (see Text S1 in the Supporting Material).

In the WLC model, harmonic potentials are used in representing the elasticities of DNABates et al. (2005); Marko and Siggia (1995b); Klenin et al. (1998). The twist energy density is defined as

where $k_{B}$ is the Boltzmann constant, $T$ is the room temperature, $l_{t}$ is the torsional persistence length, i.e., $l_{t}=75nm$, $\Delta\theta$ is the twist on the segment $\Delta u$. We introduce an internal torque $T_{q}^{in}$ as a measure of torsional stress on DNA, $T_{q}^{in}\equiv k_{B}Tl_{t}\frac{\Delta\theta}{\Delta u}$ (Torque is proportional to twist). Based on the elastic properties of DNA, the torque transport along DNA can be written as (see Text S2),

where drag coefficient for a rotational cylinder $\zeta_{R}=4\pi R^{2}\eta$, and the radius of DNA $R=2nm$ and the viscosity $\eta=0.001kg/m/s)$. A solution to Eq A2 can be seen in Text S2. The characteristic time for torque equilibrium time can be reached at $2.5\frac{\zeta_{R}}{l_{t}k_{B}T}L^{2}$(about 0.1$ms$ for $3000bp$) (Text S2), independent of stretching force, torque and linking number. For example, a twist pulse starting at one end of the stretched DNA or S phase will spread over $10^{3}$ $bp$ within $10^{-2}ms$ and $10^{4}$ $bp$ within $1ms$. Compared with the dynamics of plectoneme of $\sim 10^{3}$ $bp$ that happens at the time scale from $ms$ to $s$, e.g. diffusion and hopping, the torque transport within the S phase is fast.

Fig.1 shows the scheme for simulating the torque transport on a stretched DNA using the dWLC method. Under a constant stretching force $f=0.4$ $pN$, DNA is twisted by an external torque $T_{q}^{ex}=9pNnm$ starting at $t=0$ ($T_{q}^{ex}=0$ at $t<0$). After that, the torque is kept constant and the torsional response of DNA is monitored.

The measurements of the mean internal torque of DNA at lengths of 1500 $bp$ and 3000 $bp$ are shown in Fig. 1. After being initially twisted, the mean internal torque, $\bar{T}_{q}^{in}$, reaches to equilibrium (9 $pNnm$) at about $0.03ms$ for 1500 $bp$ and about $0.1$ $ms$ for 3000 $bp$, consistent with the torque equilibrium time $2.5\frac{\zeta_{R}}{l_{t}k_{B}T}L^{2}$.

Supercoiling can be quantified by the excess linking number $\Delta Lk$ or the linking number density $\sigma$ (see Text S3). The free energy density of a stretched DNA is, $\mathcal{F}_{S}=\frac{1}{2}c_{S}(\sigma_{S}^{eq})^{2}+\mathcal{F}_{0}$(Text S4), including the contributions from torsion and force stretchingMarko (1998, 2007), where $c_{S}$ stands for the torsional stiffness, which is a function of $f$, $l_{t}$ and $l_{b}$ ($l_{b}=50nm$, the DNA persistence length), $\sigma_{S}^{eq}$ is the equilibrium linking number density of the stretched phase, and$\mathcal{F}_{0}$ is the free energy density without torsion, representing the stretching of the WLCMarko and Siggia (1995a). At equilibrium, the torque is defined as $\frac{1}{\Omega_{0}}\frac{\partial\mathcal{F}_{S}}{\partial\sigma_{S}^{eq}}$Marko (2007), where the constant $\Omega_{0}=1.85nm^{-1}$ is the rotation angle of the DNA backbone per unit lengthBates et al. (2005); Marko (2007). Thus, we have

We accordingly measure the excess linking number density of the stretched DNA as shown in Fig. 1. $\sigma_{S}$ reaches to the equilibrium value at 0.03$ms$ for 1500 $bp$ and 0.1 $ms$ for 3000 $bp$.

A plectoneme can be considered as a pair of rods anti-parallel interwound around each other under an external torqueRybenkov et al. (1997); Marko (2007). Thus, an empirical form of the energy of the plectonemeKlenin et al. (1991); Rybenkov et al. (1997); Marko (2007) can be utilized,

where $c_{P}\equiv l_{p}k_{B}T\Omega_{0}^{2}$ is the torsional stiffness of plectoneme, and $l_{p}$ is the torsional persistence length of the plectonemic helix, depending on the ionic strength and often taking value of 21 to 27 $nm$ Rybenkov et al. (1997), and $\sigma_{P}$ is the excess linking number density of the plectonemes. The internal toque is defined similarly as that of the stretched DNA, i.e., $T_{q}^{in}=k_{B}Tl_{p}\sigma_{P}$. Nevertheless, torque transport on plectoneme can not be described as the twist angle $\theta$ propagation on DNA. That is because that the change of twist of the plectomeical helix also influences the writhe, and in turn the writhe change feeds back to the twist propagation. We introduce a plectonemic angle $\Theta$ (see Fig. S1), similar to $\theta$, as $T_{q}^{in}=k_{B}Tl_{p}\frac{\Delta\Theta}{\Delta u}$. Thus we consider the torque transport on plectoneme as the propagation of the plectonemical angle, i.e.,

where $\zeta_{P}\approx 2\pi R_{plect}^{2}\eta$ is the drag coefficient due to motions of the two parallel rods and $R_{plect}$ the radius of plectonemesNeukirch and Marko (2011); van Loenhout et al. (2012), about 2.5 nm. The time for torque equilibrium on plectoneme with contour length $L$ can be calculated as $2.5\frac{\zeta_{P}}{l_{p}k_{B}T}\frac{L^{2}}{4}$, about $2.6$ times comparing to that on the stretched DNA. i.e, about 0.26 $ms$ for torque equilibrium for 3000 $bp$.

Fig. 1 illustrates the setup of torque transport on a pure plectoneme, which is originally maintained by a pre-existing external torque $T_{q}^{ex}=8pNnm$($t<0$). Then the external torque jumps to $T_{q}^{ex}=10pNnm$ at $t\geqslant 0$, and the ensuing torsional response of the plectoneme is monitored.

The measurements of the internal torque of plectonemes at lengths of 1500 $bp$ and 3000 $bp$ are shown in Fig. 1. Following the torque jumping at $t=0$, the mean internal torques, $\bar{T}_{q}^{in}$ rises from $8pNnm$ to $10pNnm$ at about 0.1 and 0.3 $ms$ for 1500 $bp$ and 3000$bp$, respectively, as estimated above.

In equilibrium, Eq. A4 becomes $\mathcal{F}_{P}=\frac{1}{2}c_{P}(\sigma_{P}^{eq})^{2}$. The equilibrium torque is defined as $\frac{1}{\Omega_{0}}\frac{\partial\mathcal{F}_{P}}{\partial\sigma_{P}^{eq}}$. Thus

We also numerically measure the excess linking number density of the plectoneme as shown in Fig. 1. As Eq A6 suggests, $\sigma_{P}^{eq}=\frac{T_{q}^{ex}\Omega_{0}}{c_{P}}=0.076$ at $t<0$, where $T_{q}^{ex}=8pNnm$. Then $\sigma_{P}$ follows the external torque jumping. Subsequently, $\sigma_{P}$ reaches a new equilibrium defined by $\frac{T_{q}^{ex}\Omega_{0}}{c_{P}}=0.091$ at about $0.1$ ms for 1500 $bp$ and 0.3 ms for 3000 $bp$.

The torque equilibrium (Eqs A3 and A6) can be applied to the coexistence state of S and P phases. The interphase torque equilibrium means the equilibrium at the phase-boundaries, i.e., $c_{P}\sigma_{P}^{eq}=c_{S}\sigma_{S}^{eq}$. Next, we show that the phase boundaries can be chosen as the slow observables to quantify the plectoneme dynamics.

As the torque equilibrium is reached, i.e., $\Omega_{0}\overline{T}_{q}^{in}=c_{P}\sigma_{P}^{eq}$. The torque equilibrium at the boundaries means $c_{P}\sigma_{P}^{eq}=c_{S}\sigma_{S}^{eq}$. Accordingly, every plectenome carries the equal $\sigma_{P}^{eq}$. The free energy density of each plectoneme is therefore represented as $\mathcal{F}_{P}(\mathbf{X})=\frac{1}{2}c_{P}(\sigma_{P}^{eq})^{2}$Klenin et al. (1991); Rybenkov et al. (1997).

The translational displacement of the $\alpha$-th section of the S phase is subjected to the viscous drag along it. At $u<X_{\alpha}^{r}$ and $>X_{\alpha}^{l}$ , the tension imposed on DNA is $f_{u}$, and at the stretching end, $f|_{u=L}=f$. As the torque equilibrium is reached, the local free energy density (or force) is represented by $\mathcal{F}_{S}(u,\mathbf{X})=\frac{1}{2}c_{S}(\sigma_{S}^{eq})^{2}+\mathcal{F}_{0}$Marko (1998, 2007).

The free energy associated with ${\mathbf{X}}$ can be written as the contributions from the S and P phases

The free energy gradient along $\mathbf{X}$ serves as a driving force, which drives the transformation between S and P phases. The variable $\mathbf{X}$ thus follows the Langevin equations $\gamma_{1}\mathbf{\dot{X}}=-\nabla_{\mathbf{X}}\Phi_{0}(\mathbf{X})+\sqrt{\frac{2k_{B}T\gamma_{1}}{dt}}\mathbf{\dot{W}}$, where $\gamma_{1}$ denotes the drag coefficent for the growth or shrinkage of plectonemes, and $\dot{W}$ is the white noise. Specifically, the left/right boundaries of the $\alpha$-plectoneme follow

where $\gamma_{\alpha}$ is drag coefficient for the growth or shrinkage of the $\alpha$-th plectoneme (see Methods).

It should be noted that Eq B2 should be completed with two constraints: the first one is the topology constraint ($\Delta Lk=\Delta Lk_{S}+\Delta Lk_{P}.$), and the second one is in-extensibility of DNA model ($L=L_{S}+L_{P}$)(see Methods).

The equilibrium point, i.e., $\nabla_{X_{\alpha}^{l/r}}\Phi_{0}(\mathbf{X})=0$ leads to the coexistence torque $T_{co}=\frac{1}{\Omega_{0}}(\frac{2c_{P}(f-\sqrt{\frac{k_{B}Tf}{l_{b}}})}{1-c_{P}/c_{S}})^{\frac{1}{2}}$Marko (2007). Actually, the dynamics of the slow variable $\mathbf{X}$ accompanies with a time-dependent torque, $T_{q}^{in}(t)\equiv\frac{c_{P}\sigma_{P}^{eq}(t)}{\Omega_{o}}$ (Fig. 2). The internal torque of a DNA at length of 4500 $bp$ under $f=0.5pN$ fluctuates around the coexistence torque, $T_{co}$. The two-phase dynamic model can be regarded as a smoothed version of the dWLC method. Although the fluctuations faster than the torque equilibrium have been averaged out, the mean and the variance of the torque are still consistent with that obtained from the dWLC method(inset in Fig. 2).

Supercoil with a given linking number is analogous to the classical gas (van der Waals) confined inside a fixed volume (NVT ensemble)Mattis and Swendsen (2008). A DNA twisted under a constant torque $T_{q}^{ex}$, however, is analogous to the classical gas under a fixed pressure (isobaric ensemble). Its free energy densities can be obtained through the Legendre’s transformation $\mathcal{F}_{S/P}(T_{q}^{ex})=\mathcal{F}_{S/P}-\Omega_{0}T_{q}^{ex}\sigma^{eq}$. The free energy densities of the plectonemes and the stretched phase are therefore $\mathcal{F}_{P}=-\frac{1}{2}\frac{(\Omega_{0}T_{q}^{ex})^{2}}{c_{P}}$ and $\mathcal{F}_{S}=-\frac{1}{2}\frac{(\Omega_{0}T_{q}^{ex})^{2}}{c_{S}}-f_{u}+\sqrt{\frac{k_{B}Tf_{u}}{l_{b}}}$, respectively. Inserting them into Eq B1 and Eq B2 we obtained the dynamic model of extended DNA supercoiling under constant torque.

Generally speaking, the stretching force, $f$ in Eq B1 or B2 is not necessarily constant. The force can be time-dependent as long as it dose not break the interphase torque equilibrium. For example, one can impose an instantaneous force. Similarly, $\Delta Lk$ can be time-dependent as well if it does not break the interphase torque equilibrium.

An extended DNA undergoes linking number accumulation when one end of the DNA is rotated while the other end is kept fixed, as shown in Fig.3. The stretching force $f=0.5pN$, the rotation angle $\Omega=\omega t$ and the angular velocity $\omega=20\pi/s$, equivalent to $\Delta Lk=10$ per second (i.e., comparable to polymerase enzyme unwinding and synthesizing at  100 nucleotide or nt per second)Thomen et al. (2005); Jing et al. (2018). At first, the extended DNA is twisted, and then it undergoes a buckling transition (discontinuities on torque and extension curves) followed by plectoneme formationForth et al. (2008).

We performed BD simulations of such a supercoil accumulation based on the dWLC method. In Fig. 3, the relative extension $z/L$ almost remains unchanged before 0.6 s while the torque gradually builds up, suggesting that the extended DNA is twisted. The buckling at about 0.6 s marks the starting-point of the phase-coexistence. The buckling transition indicates a competition between pre-plectonemic loops and stretched phaseMarko and Neukirch (2012); Mielke et al. (2004). After 0.6 s, the torque drops and converges as torsional stress on DNA is relieved, indicating the twist is transformed into writhe. These discontinuities have previously been observed in single-molecule experimentForth et al. (2008). The coexistence torque $T_{q}^{co}$ remains unchanged when more linking number is injected. Note that multiple plectonemes can co-exist in the process. The representative snapshots are also shown in Fig. 3. The unchanged torque curves or the linearly decreased extension after buckling transition suggest that the process remains quasi-static.

It should be noted here that we obtained the curves in Fig. 3 from four BD simulation trajectories using the dWLC method, and each of the trajectories requires several weeks of CPU time. The calculation of the electrostatics interaction between the charged points on segments of the dWLC model indeed dominates the computational cost (see Text S6). According to our analyses, the computational cost of the dWLC in the plectonemic phase is still significant and grows linearly with the contour-length. For example, in Fig. 1, simulating a stretched DNA of 3000 $bp$ in 10 $ms$ costs 3 hours while simulating a pure plectoneme of the same contour length costs one day(i.e., the CPU time). The two-phase dynamic model, however, enables us to reduce the computational cost to trivial by dealing only with the essential variables with slow dynamics. In Fig. 3, the torque and extension curves generated from the two-phase dynamic model (using tens of seconds of computation time) well reproduce those in Fig.3(b) from the dWLC method and also consistent with previous experimental and computational studiesForth et al. (2008); Ott et al. (2020).

We also performed simulations of supercoiling generation under a constant torque. Fig.3 illustrates an extended supercoil under a constant torque, $T_{q}^{ex}=9pNnm$ (comparable to 5$\sim$11 $pNnm$ torque exerted by polymerasesMa et al. (2013)), and a constant stretching force $f=0.4pN$. Since the external torque is larger than $T_{co}$ (i.e., $6.5pNnm$ that maintains the coexistence state under 0.4 $pN$), plectonemes are created and the ultimate equilibrium state is dominated by plectonemes. Consequently, plectonemic coils are building up until all stretched DNA is interwound to plectonemical coils.

Using the dWLC method, we can capture the dynamics of the plectoneme growth. Fig. 3 shows the measurements of the relative extensions $z/L$ of DNAs at 1500 $bp$ (black solid squares), 3000 $bp$, 6000 $bp$ and 9000 $bp$. Indeed, the shrinking extensions indicate the growth of plectonemes. The fitting curves are generated using $\frac{z}{L}-\frac{z_{0}}{L}\propto 1-\sqrt{1+\lambda t}$, where $\lambda$ is a fitting parameter. In the two-phase dynamic model, all but the slithering drag coefficients are given based on the hydrodynamics (see Methods). The slithering drag $\mu_{slt}$ is to be determined in the two-phase model by calibrating results with the dWLC method. We performed the simulations using the two-phase dynamic model shown in Fig. 3. The slithering drag $\mu_{slt}$ was then determined via fitting parameter $\lambda$ Text S7). Consequently, we can then reproduce the plectoneme dynamics.

The inter-phase boundaries in the two-phase dynamic model can be used to directly monitor the plectoneme dynamics. Fig. 4 shows a kymograph, i.e., the spatial vs time evolution process of the plectonemes of a 21-kbp (7$\mu m$) DNA (extended under a force $f=0.8pN$ and with $\Delta Lk=70$), illustrating the coexistence of multiple plectonemes, and the plectoneme nucleating and vanishing.

The simulations of a plectoneme with $1/4$ of the DNA contour length under different forces were performed. To estimate the diffusion constant, we measured the mean square displacement (MSD) $\overline{\Delta x^{2}(t)}$ of the plectoneme. As shown in Fig. 4, we obtained the diffusion constants via $D\equiv\frac{\overline{\Delta x^{2}(t)}}{2t}$ under different stretching forces. For consistency check, we also compared the results with the Einstein relation, i.e., $D=\frac{k_{B}T}{\zeta_{diff}}$ where $\zeta_{diff}$ as a sum of the viscous drags for the transverse displacement of the plectoneme (see Methods). The diffusion constant of the plectoneme is about $10^{-1}\mu m^{2}/s$.

To compare the dynamics of the plectonemes at different stretching forces and ionic strengths, simulations of the P phase with $1/4$ of the DNA contour length were performed in the conditions where the stretching forces were $0.4$, $0.8$, $1.6$ and $3.2$ $pN$. Large stretching force $f$ results in more compact plectonemes, i.e., reduces the radius of the plectonemic helix(Fig 4). Indeed, in the coexistence state of DNA supercoiling, the equilibrium excess linking number density of the P phase is $\sigma_{P}^{eq}=\frac{1}{c_{P}}(\frac{-2c_{P}(f-\sqrt{\frac{k_{B}Tf}{l_{b}}})}{1-c_{P}/c_{S}})^{\frac{1}{2}}\sim\sqrt{f/l_{p}}$Marko (2007), which suggests that $\sigma_{P}^{eq}$ increases with $f$.

At high salt concentration, the screening of the Coulomb repulsion potential reduces the effective repulsive DNA diameter, leading to the shorter torsional persistence length of the plectoneme $l_{p}$ and smaller radius of the plectonemic helixRybenkov et al. (1997). Meanwhile, based on the coexistence state relation $\sigma_{P}^{eq}\propto\sqrt{f/l_{p}}$, lowering $l_{p}$ increases the $\sigma_{P}$, i.e., making the plectonemes more compact (Fig. 4). Accordingly, we adjusted the torsional persistence length of plectoneme $l_{p}$ to represent the impacts from the ionic strengthRybenkov et al. (1997), specifically, $l_{p}=25$, $23$ and $21$ $nm$ to represent the low, regular and high salt concentrations, respectively.

According to Kramers’ reaction-rate theory, the nucleation rate is $k_{nuc}\propto exp(-\beta\epsilon_{nuc})$ Daniels and Sethna (2011), where the energy barrier $\epsilon_{nuc}\propto\sqrt{f}$, which incorporates the bending energy of the loop and the loss of the elastic energy stored originally in the stretched DNA(Fig. 4)Marko and Neukirch (2012). Therefore, the large forces can lower the nucleation rate. Again, according to the relation $\sigma_{P}^{eq}\propto\sqrt{f/l_{p}}$, large force and small $l_{p}$ (high salt concentration) have a similar effect on the nucleation rate. Indeed, the nucleation rates obtained from simulations (Fig. 4) suggest that large stretching forces $f$ and the high salt concentrations (small $l_{p}$) can suppress the nucleation events. We derived the nucleation relation (Text S8)

where $2C$ and $l_{p}^{\star}$ are fitting parameters (dashed curves in Fig. 4). As $l_{p}\rightarrow l_{p}^{\star}$ (salt concentration decreases), the nucleation rate increases. That is because low salt concentration increases the radius of the plectonemic helices i.e.,$R_{plect}\rightarrow R_{loop}$ (bending energy of DNA in plectonemes increases), then lower the energy barrier from plectonemic helix to loop (reversed buckling transition, red arrow in Fig. 4). By fitting, we obtained $l_{p}^{\star}=27nm$.

The results above indicate that the plectonemes frequently vanish and re-emerge (via nucleation), which lead to the hopping. The surviving time of the individual plectoneme from nucleation to vanishing is defined as plectoneme life-time. Fig. 4 shows that the mean life-time $\tau_{plec}$ depends on stretching force and salt concentration. The large $f$ or small $l_{p}$ increases the plectoneme life-time. As discussed above, large $f$ or small $l_{p}$ lowers the plectonemes nucleation and also leads to more compact plectonemes. The resulting compactness stabilizes existing plectonemes, which means long mean life-time of the individual plectonemes. By obtaining the balance between the plectoneme nucleation and vanishing, we can derive a mean life time of the individual plectoneme (see Text S8)

where $C$ and $l_{p}^{\star}$ are the same as those in Eq D1.

An example of hopping $over10kbp$ within about $0.1s$ is shown in Fig.4. Plectoneme hopping is a special case of the DNA conformational rearrangementvan Loenhout et al. (2012). As addressed above, plectoneme diffusion happens comparatively slowly, only hundreds of bps per second ($\sim\sqrt{2Dt}$). Hopping , however, is a rapid process, potentially, facilitating fast site-specific recombination and enhancer activated expressionSheinin and Wang (2012). Since plectoneme hopping is related to nucleation, the large forces inhibiting the nucleation thus repress the hopping.

We also found that plectoneme hopping is associated with significant energetic changes between the two phases. We identify the nucleation as the starting point of a hopping event. Then we collected numerous simulations (5000) that contain the hopping events and aligned the starting points of the hopping events together for averaging the energetics. The energy change of the S phase accumulates up to about $20k_{B}T$ before the nucleation, while that of the P phase falls down to about $20k_{B}T$. These significant energetic changes facilitate the abrupt nucleating-vanishing processes (Fig. 4). After the formation of the new plectoneme by nucleation, the energetic changes rapidly converge to $0$. The energetic offsets between the two phases thus makes the total energy change smaller as thermal fluctuations, i.e.,$1k_{B}T$(inset of Fig. 4). It should be noted that there is no persistent energetic injection into the system. Thus the energy change must be compensated by the complementary reversed processes (one of the multiple plectonemes vanishes, corresponding to reversed buckling transition shown in Fig. 4).