Transcription-driven DNA Supercoiling: Non-Equilibrium Dynamics and Action-at-a-distance

TL;DR

This study uses simulations to reveal how transcription causes asymmetric DNA supercoiling and influences DNA-protein interactions at a distance, challenging traditional symmetric models.

Contribution

It introduces a detailed simulation framework showing asymmetric supercoiling domains and action-at-a-distance effects of transcription on DNA structure.

Findings

Supercoiling domains are asymmetric, with plectonemes forming far from RNAP.

Positive supercoils destabilize nucleosomes before RNAP reaches them.

Twist diffusion is fast, while writhe relaxation is slow and multi-step.

Abstract

We study the effect of transcription on the kinetics of DNA supercoiling in 3D by means of Brownian dynamics simulations of a single nucleotide resolution coarse-grained model for double stranded DNA. By accounting for the action of a transcribing RNA polymerase (RNAP), we characterise the geometry and non equilibrium dynamics of the twin supercoiling domains forming on each side of the RNAP. Textbook pictures depict such domains as symmetric, with plectonemes (writhed DNA) appearing close to the RNAP. On the contrary, we find that the twist generated by transcription results in asymmetric domains, with plectonemes formed far from the RNAP. We show that this translates into an "action-at-a-distance" on DNA-binding proteins: for instance, positive supercoils downstream of an elongating RNAP destabilise nucleosomes long before the transcriptional machinery reaches the histone octamer. To understand these observations we use our framework to quantitatively analyse the relaxation dynamics of supercoiled DNA. We find a striking separation of timescales between twist diffusion, which is a simple and fast process, and writhe relaxation, which is slow and entails multiple steps.

Figures

b in the main text for example). We then remove both the topological barrier and the RNAP complex, and we analyse the evolution of this new system to its relaxed state. Under these conditions, the two opposite sign plectonemes will relax and eventually disappear. One way of studying this relaxation process is by defining an effective local writhe field $\chi(n,t)$ at position $n$ along the DNA and time $t$, that can be tracked directly from our Brownian dynamics (BD) simulations. To define $\chi(n,t)$ we consider a segment of the DNA as a “sliding window” centred at bp $n$. The length of this

c and Fig. [S9](#S13.F9)) using a phenomenological ODE,

Equations

$$U_{\mathrm{LJ}}(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6+\frac{1}{4}\right],$$

$$U_{\mathrm{bb}}(r)=\{\begin{array}{ll}-\frac{{ϵ}_{bb}{R}_0^2}{2}\text{ln}\left[1-{\left(\frac{r}{R_0}\right)}^2\right]& \text{if }r<R_0\\ \infty ,& \text{if }r\ge R_0,\end{array}$$

$$U_{\mathrm{morse}}(r)={ϵ}_m[1-e^{-\lambda (r-r_{0\text{,stack}})}{]}^2,$$

$$U_{\mathrm{planarity}}(\alpha )=\frac{{ϵ}_p}{2}(\alpha -{\alpha }_0{)}^2,$$

$$U_{\mathrm{dihedral}}(\varphi )={ϵ}_d[1+\text{cos}(\varphi -d)],$$

$$U_{hb}(\mathrm{r})=\frac{{ϵ}_{\mathrm{HB}}}{({\mathrm{r}}_{0\text{,hb}}-{\mathrm{r}}_{\text{c,hb}}{)}^2}[(\mathrm{r}-{\mathrm{r}}_{0\text{,hb}}{)}^2-({\mathrm{r}}_{\mathrm{c}\text{,hb}}-{\mathrm{r}}_{0\text{,hb}}{)}^2],$$

$$U_{\mathrm{bending}}={ϵ}_b[1+\text{cos}(\theta )].$$

$$m\frac{d^2\mathbf{r}}{d{t}^2}=-\xi \frac{d\mathbf{r}}{dt}-\nabla \mathcal{U}+\sqrt{2k_{B}T\xi }\eta (t),$$

$$⟨{\eta }_{\alpha }(t){\eta }_{\beta }(s)⟩={\delta }_{\alpha \beta }\delta (s-t)$$

$${\tau }_{\mathrm{in}\mathrm{rot}\mathrm{solid}\mathrm{sphere}}=\frac{3}{10}{\tau }_{\mathrm{in}}.$$

$$\sin \theta (n)=\mathbf{t}(n)\cdot [{\mathbf{m}}^{\ast }(n)\times \mathbf{m}(n)].$$

$$\mathrm{Tw}=\frac{1}{2\pi }\sum _n\theta (n).$$

$$\mathrm{Wr}=\frac{1}{4\pi }{\int }_C{\int }_C\frac{{\mathbf{r}}_2-{\mathbf{r}}_1}{|{\mathbf{r}}_2-{\mathbf{r}}_1|}d{\mathbf{r}}_2\times d{\mathbf{r}}_1,$$

$$\mathrm{Wr}=2\sum _{i=2}^N\sum _{j<i}\frac{{\Omega }_{ij}}{4\pi },$$

$${\mathbf{n}}_i$$

$${\mathbf{n}}_{i+1}$$

$${\mathbf{n}}_j$$

$${\mathbf{n}}_{j+1}$$

$${\Omega }_{ij}^{\ast }$$

$$\arcsin ({\mathbf{n}}_j\cdot {\mathbf{n}}_{j+1})+\arcsin ({\mathbf{n}}_{j+1}\cdot {\mathbf{n}}_i).$$

$$\frac{{\Omega }_{ij}}{4\pi }=\frac{{\Omega }_{ij}^{\ast }}{4\pi }\text{sign}[({\mathbf{t}}_{j,j+1}\times {\mathbf{t}}_{i,i+1})\cdot {\mathbf{t}}_{i,j}],$$

$$\frac{\partial \theta (x,t)}{\partial t}=D_{\mathrm{Tw}}\frac{{\partial }^2\theta (x,t)}{\partial x^2},$$

$$\theta (x,0)={\theta }_0+\Delta \theta \text{H}(x-x_i)\text{H}(x_f-x).$$

$$\theta (x,t)=\{\begin{array}{l}\frac{{\theta }_0}{2}\left[2+\text{erf}\left(\frac{x-x_i}{2\sqrt{Dt}}\right)-\text{erf}\left(\frac{x-x_f}{2\sqrt{Dt}}\right)\right]\text{, for }\Delta \theta >0\\ \\ \frac{{\theta }_0}{2}\left[2+\text{erf}\left(\frac{x-x_f}{2\sqrt{Dt}}\right)-\text{erf}\left(\frac{x-x_i}{2\sqrt{Dt}}\right)\right]\text{, for }\Delta \theta <0,\end{array}$$

$$\frac{E_{\mathrm{bend}}}{k_{B}T}=\frac{l_p}{2l}\sum _{i=0}^N|{\mathbf{t}}_{i+1}-{\mathbf{t}}_i{|}^2,$$

$$\frac{\partial w(t)}{\partial t}=-D\left(w(t)\right)w(t),$$

$$D(w(t))=d_1+(d_2-d_1)\frac{1+\tanh [b(w(t)-w_c)]}{2}.$$

Full text

Keywords— DNA topology $|$ Supercoiling $|$ Transcription $|$ Non-Equilibrium Physics

The double stranded nature of DNA endows it with the properties of an elastic rod with both bending and twisting rigidities [1, 2]. Mechanical manipulation which over or under winds the double helix leads to torsional strain that can be relieved by the DNA writhing onto itself [2]. This phenomenon is a consequence of a topological conservation law: if the DNA is in a closed loop, or its ends are fixed, the number of times the two strands wind around each other – the linking number ($\mathrm{Lk}$) – is a topological invariant. The White-Fuller-Calugarenau (WFC) theorem [3, 4, 5] states that the linking number of a DNA helix can be written as the sum of two contributions, its “twist” ($\mathrm{Tw}$) and “writhe” (${\mathrm{W}r}$), $\mathrm{Lk}=\mathrm{Tw}+\mathrm{Wr}$. The twist can be thought of as the number of times the vector joining the two DNA nucleotides in a base pair rotates around the backbone, whereas the writhe counts the (signed) self-crossings of the backbone [2]. In its relaxed state a DNA double helix has a unit linking number for every $\sim 10.5$ base-pairs (bps); a molecule with a linking number different from this is said to be supercoiled.

Bacterial plasmids are kept in a negatively supercoiled state (corresponding to an underwound double helix), and this is thought to facilitate transcription [6]. In eukaryotes, the modulation of DNA twist along the chromosome is thought to play a key role in gene regulation [7]. Irrespective of the organism, in vivo DNA is constantly being remodelled by proteins such as RNA polymerase (RNAP), that drive the system away from equilibrium by applying forces and torques of the order of $25$ pN [8, 9] and $11$ pN nm [10], respectively. Force and torque generation is required for transcription, and inevitably introduces supercoiling in a DNA molecule – this fact is the basis of the “twin supercoiling domain” (TSD) model [11, 12]. The TSD model predicts the formation of positive supercoils ahead of the RNAP and negative supercoils in its wake. Such dynamically generated supercoiling has been conjectured to play a regulatory role in gene expression either through the twist dependence of polymerase-DNA interactions [13, 14], or via supercoiling-mediated generation of DNA loops [6]; it may also affect chromatin architecture close to promoters [15]. Although the TSD model was proposed over 30 years ago, its consequences in vivo are still far from understood.

In this work we use 3-D Brownian dynamics (BD) simulations of a single nucleotide resolution model for DNA [16, 17] to study the transcriptionally-induced generation of TSDs. Our simulations go beyond the 1-D description of RNAP-driven supercoiling used in recent models [13, 14, 18] as well as earlier studies using a twistable worm-like chain model for DNA [19]. Contrary to typical textbook diagrams which show twist or writhe being generated close to the RNAP on both sides, we discover that writhe nucleates into plectonemes that appear at a considerable distance from the polymerase: this allows for an “action-at-a-distance” phenomenon, where transcription at one point on the DNA can affect its dynamics and the binding of proteins – such as histones or other RNAPs – elsewhere. Additionally, we study twist and writhe relaxation in the absence of RNAP and find they differ vastly, both quantitatively and qualitatively: twist diffuses away rapidly, whereas writhe relaxation is much slower and entails at least two distinct timescales. Our results challenge the TSD paradigm and motivate further effort to dissect the dynamics of supercoiling in single-molecule experiments.

Results

The Model

We use a recently developed single nucleotide resolution model for double-stranded (ds) DNA that fully captures its double helical structure and elastic properties [16, 17], and simulate a closed DNA loop. Unless otherwise stated, we consider an $L=1000$ base-pair (bp) long loop in which the the linking number is equal to the “relaxed” value $\mathrm{Lk}=\mathrm{Lk}_{0}=L/10.5$.

An RNAP is modelled as a rigid body consisting of a ring which encircles the double helix and a “crossbar” which passes between the two DNA strands and whose ends are anchored to the ring (see Fig. 1a). The dsDNA segment passing through the RNAP experiences an active force, $F_{p}$, directed perpendicular to the plane of the ring, which drives the relative motion of RNAP and DNA. The steric interaction between the crossbar and the DNA beads forces the opening of the double strand and leads to positive twisting (overwinding) in front of the RNAP and negative twisting (underwinding) behind. Our model correctly captures the mechanical action of an elongating RNAP which must unwind a section of the helix in order to “read” the nucleotides [20, 1]. We quote forces in units of $\epsilon/\sigma$, where $\epsilon$ and $\sigma$ are the simulation energy and length units respectively; full details of the simulation units, and the DNA and RNAP models are given in SI Appendix.

For simplicity, here we focus on a geometry where the RNAP is fixed in 3-D space; this could either mimic an in vitro set-up where the RNAP is fixed in place, or capture the fact that in vivo RNAPs are associated with a large elongation complex which experiences a larger rotational drag than the DNA [21]. [Similar results are obtained for a moving RNAP, SI Appendix and Fig. S3).]

RNAP Elongation Creates Asymmetric TSDs

Within the crowded environment of the cell, binding of large protein complexes or the formation of loops can restrict DNA rotation at specific sites, hindering the transmission of twist [20]. To mimic this constraint we introduce a “topological barrier” (TB), a section of DNA which is not able to rotate and thus effectively acts as a reflecting barrier with respect to twist transport (Fig. 1a). In this context, we investigate the effect of twist generation by examining different values of $F_{p}$ for a DNA molecule which was initially topologically relaxed (Fig. 1b, left).

Once a sufficiently large amount of supercoiling has been introduced by the RNAP, we observe the emergence of plectonemic structures on each side. Unexpectedly, we discover that plectonemes can appear at large distances from the RNAP complex, and that each can store a different amount of writhe. For instance, to the right in Fig. 1b we show a configuration obtained after the RNAP has travelled $3.4$ turns of the helix: plectonemes with $\mathrm{Wr}=+1$ and $\mathrm{Wr}=-2$ have formed ahead of and behind the RNAP respectively. This is contrary to the usual TSD picture, which would suggest equal and opposite writhing on each side close to the RNAP [11]. Even though the writhe within the plectonemes is unbalanced, the total linking number must still be conserved to satisfy the WFC theorem, i.e. $\mathrm{Lk}(t)=\mathrm{Lk}_{0}=\mathrm{Tw}(t)+\mathrm{Wr}(t)$ [2]. We verify that this is the case by computing the global twist and writhe as a function of time (Fig. S4; see SI Appendix for details on how these quantities are computed). This indicates that not all of the writhe is stored in the plectonemes, but part of it is delocalised over the whole polymer conformation.

Plectonemes Form At a Distance From the RNAP

To understand the formation of plectonemes more quantitatively, we developed a strategy to monitor both their position and length as they are generated (see Fig. S5). Briefly, at each time step, we compute the matrix of contacts between each DNA bp and generate a simulated contact map [22]. Two non-adjacent bps come into close proximity only at the crossing points within plectonemes and we define the length of a plectoneme, $l$, as the contour length between the segments forming the outermost crossing. We also measure the distance $d$ between the RNAP and the base of the first plectoneme. Further details are given in SI Appendix and Fig. S5.

First, we find that the larger the force the smaller the plectoneme (Fig. 1d); specifically $l\sim F_{p}^{-\alpha}$ with $\alpha\simeq 0.9$–$1$. This scaling differs from the $\sim F_{p}^{-1/2}$ observed for plectoneme formation in stretched DNA at equilibrium [23], and is compatible with a Pincus scaling for blob size $\xi\sim k_{B}T/F_{p}$. [We expect the tension to be non-uniform in our set-up, as the force is applied locally at the polymerase location only; however the local tension should still scale with $F_{p}$ [24].]

Second, we find that depending on $F_{p}$ there are two regimes for the distance $d$ between the RNAP and the plectoneme (Fig.1e). For large forces ($F_{p}\geq 30~{}\epsilon/\sigma$), $d$ depends on the sign of the supercoiling: for $\mathrm{Wr}>0$, i.e. in front of the RNAP, $d$ increases with the force, whereas for $\mathrm{Wr}<0$, i.e. behind the RNAP, it decreases with $F_{p}$. We attribute this to the tension experienced by the DNA in front of the RNAP as it is “reeled in”: the backbone is straightened and writhing is suppressed (Fig. 1e, rightmost inset). At the same time, there is an accumulation of DNA under compression behind the RNAP, and this can more readily form a plectoneme.

For small forces ($F_{p}<30~{}\epsilon/\sigma$), we observe a markedly different behaviour: the RNAP does not exert enough force to immediately break the attraction between the DNA strands, and it moves much more slowly (Fig. S6). This is because RNAP progression now depends on thermal fluctuations which favour DNA opening. The switch between the low and high force regime is reminiscent of the transition between stick-slip and sliding motion in the Prandtl-Tomlinson model for dry friction [25]. In the low force regime, the DNA is not under significant tension (or compression), and plectonemes form symmetrically. Plectoneme formation is hindered close to the RNAP and the TB due to steric effects and reduced DNA mobility, and as a result we observe DNA writhing close to the halfway point between the two (Fig. 1e).

Finally, we note that in about $80\%$ of the simulations, and irrespective of $F_{p}$, the negatively supercoiled plectoneme is the first to form. This finding is in agreement with the observation in Ref. [26] that negatively supercoiled DNA is more difficult to twist (but easier to bend) compared to positively supercoiled DNA. Furthermore it implies that our model correctly captures the asymmetric torsional rigidity of dsDNA [27].

Transcription Unwraps Nucleosomes At a Distance

We reasoned that the “action-at-a-distance” just described might mechanically affect the binding of a protein on DNA. To test this hypothesis, we consider the basic building block of eukaryotic chromatin, the nucleosome, and model a histone octamer as a spherical bead of diameter $10$ nm, with $20$ “sticky patches” tracing a left-handed helical path on its surface covering exactly $1.7$ turns [28]. By introducing a short range attraction between the patches and the DNA beads, modelling screened electrostatic interactions, we can readily simulate the self-assembly of a nucleosome (Fig. 2a).

We initialise the system as a torsionally relaxed loop of DNA, containing an initially inactive RNAP ($F_{p}=0$) and an octamer, positioned at opposite sides of the loop. We observed that as the DNA wraps around the octamer, the topology of the nucleosome is in agreement with the “linking number paradox” [2, 29]: while the DNA completes $1.7$ turns around the histone core, the wrapped section is also slightly over-twisted (see bar plot in Fig. 2a), resulting in a net linking number storage of about $-1$ per nucleosome [2].

Once the nucleosome has formed and the DNA is in an equilibrium configuration (Fig. 2a), we introduce a TB as before, and activate the RNAP. The TB isolates the nucleosome from supercoils coming from one direction, and we can dissect the role of negative and positive supercoiling by considering an RNAP moving towards or away from the nucleosome. In Figure 2b (right-hand plot) we plot the averaged local twist $\langle\theta(x)-\theta_{0}\rangle$ at a fixed time after RNAP activation at key locations along the DNA (at the nucleosome, at the TB, and within the supercoiled domains). We use as the baseline simulation time units the Brownian time $\tau_{Br}$ as defined in SI Appendix. Remarkably, when the nucleosome is subject to RNAP-driven positive supercoiling it becomes destabilised and the DNA associated with it unwraps long before the RNAP reaches it (Fig. 2b, left). Conversely, when the nucleosome is subject to negative supercoiling there is no unravelling, and instead the nucleosome structure becomes more stable (Figs. 2c and S7).

Our simulations not only agree with previous observations of transcription-driven nucleosome eviction in vitro [30], but also suggest that the removal of obstacles in front of an advancing RNAP may take place without any direct contact, through an “action-at-a-distance” mechanism which exploits the travelling of supercoiling along the DNA.

Two Modes of Supercoiling Relaxation: Twist Diffusion

Having observed that RNAP-driven supercoiling spreads along the DNA, we now use our 3-D simulations to dissect the mechanism behind twist and writhe relaxation. First, we consider twist relaxation. We initialise a $1000$ bp dsDNA loop in a non-relaxed state with $\mathrm{Tw}-\mathrm{Tw}_{0}=3$. This is done by fixing the twist angle contained within a short ($30$ bp) segment at $\theta(n)=72\degree$, away from the preferred value of $\theta_{0}=36\degree$; we subsequently run a simulation that preserves this constraint to generate an equilibrated conformation with a locally overtwisted segment. By releasing the constraint we then study the relaxation of the twist by monitoring its local value along the molecule (see SI Appendix) and averaging this over $1000$ independent simulations (Fig. 3a).

We find that twist relaxation can be fitted by an analytical solution of the diffusion equation, $\partial_{t}\theta(x,t)=D_{\rm Tw}\nabla^{2}\theta(x,t)$, with initial condition $\theta(x,0)=2\theta_{0}=72\degree$ if $x_{i}<x<x_{f}$ and $\theta(x,0)=\theta_{0}=36\degree$ otherwise (see SI Appendix). By fitting $\theta(x,t)$ to our simulation data we extract the diffusion coefficient $D_{\rm Tw}^{+}=478\pm 9$ bp${}^{2}/\tau_{\mathrm{Br}}$. A similar but slightly smaller value of $D_{\rm Tw}^{-}=312\pm 12$ bp${}^{2}/\tau_{\mathrm{Br}}$, was obtained for an under-twisted DNA, initialized with $\theta_{n}=0\degree$ within the same segment (between points $x_{i}$ and $x_{f}$). These numerical estimates are robust for different lengths of the over/under-twisted segment.

We highlight that our results suggest that overtwist diffuses faster than undertwist. This is reasonable since the two deformations are not related by a simple symmetry and because our initial local twist deformations are large (they are outside the linear regime where the behaviour of over/under-twist should be symmetric [31, 17]).

Two Modes of Supercoiling Relaxation: Glassy Writhe Dynamics

To study the relaxation of writhe, we consider initial configurations which display stable plectonemic structures (localized writhe) and a twist which is close to the relaxed value ($\theta(n)\simeq\theta_{0}$). To achieve this, we use configurations generated in our TSD simulations (such as shown to the right in Fig. 1b); we remove both the RNAP and TB and then monitor the local writhe $\chi(n,t)$ [32, 33] (see SI Appendix) as the molecule relaxes. Typical profiles for $\chi(n,t)$ at different times $t>0$ are shown in Figure 3b: the peaks and troughs identify positively and negatively writhed plectoneme tips respectively, and we can quantify their evolution by recording the maximum and minimum values of $\chi(n,t)$ [32].

We discover that the relaxation of plectonemes occurs in two steps: an initial fast relaxation followed by a slower decay at longer times (Fig. 3c). We reason that this is due to high levels of writhe stored in the plectonemes; these carry conformational stress that is quickly released as soon as the DNA is allowed to relax (see Fig. S10). At later times, the remaining writhe becomes delocalised (as indicated by the broadening of the peaks in $\chi(n)$), which entails a lower conformational stress leading to a slower decay. [A more quantitative analysis supporting this interpretation is provided in SI Appendix, Fig. S11.] We find that the fast and slow relaxations are well fitted by two exponential decays (Fig. 3c), where the second relaxation time scale is about an order of magnitude larger than the first. The data can also be fitted by a stretched exponential decay, which would be compatible with an effective fractional diffusion equation for writhe, corresponding to sub-diffusion [34].

Our results strongly suggest that writhe and twist relaxation have profoundly different timescales, and we find that both display an intrinsic asymmetry between positive and negative supercoiling. Indeed, negatively writhed plectonemes systematically relax more slowly than those involving positive writhe. More specifically, we find that the average timescales of positive writhe relaxation are $\tau_{1}^{+}\simeq 2.9$ $10^{3}\tau_{\mathrm{Br}}$ and $\tau_{2}^{+}\simeq 29$ $10^{3}\tau_{\mathrm{Br}}$ whereas for negative writhe we measure $\tau_{1}^{-}\simeq 4.4$ $10^{3}\tau_{\mathrm{Br}}$ and $\tau_{2}^{-}\simeq 44.5$ $10^{3}\tau_{\mathrm{Br}}$ (see SI Appendix and Table S1 for details and additional results). Since there is no energy difference between a positively and negatively writhed configuration, this asymmetry points to a writhe relaxation pathway which involves an intermediate partial conversion to twist.

Together these observations suggest that supercoiling relaxation occurs through two distinct mechanisms: twist diffusion and writhe dissipation. While the former takes place on very short timescales – $D_{\rm Tw}\simeq 500$ bp${}^{2}/\tau_{\mathrm{Br}}$ – the latter is much slower. The effective diffusion coefficient for writhe can be estimated as $D_{\rm Wr}\simeq l^{2}/(\tau_{1}+\tau_{2})\sim$  bp${}^{2}/\tau_{\mathrm{Br}}$ (where $l\sim 100$ bp is the initial plectoneme size), over two orders of magnitude smaller than $D_{\rm Tw}$. This difference arises because writhe relaxation requires global conformational changes which are not necessary to dissipate twist.

Transcription without topological constraints

While it is reasonable to assume that in vivo DNA is subject to topological constraints similar to those imposed by our topological barrier, such constraints may be absent for in vitro set-ups. It is therefore relevant to ask if, in the absence of a TB, RNAP-driven supercoiling can still lead to nonequilibrium plectoneme formation, or whether it would simply travel around the DNA loop as twist and annihilate. To answer this question we measure the total unsigned writhe $\zeta(t)$ (see SI Appendix for its definition) generated by an RNAP in a DNA loop with no TB. This “order parameter” quantifies the appearance of DNA crossings, irrespectively of their sign [32, 33], and is thus non-zero only in the regime in which plectonemes form.

Our simulations show that $\zeta(t)$ reaches a steady state $\bar{\zeta}$ and, when plotted as a function of the force, $\bar{\zeta}$ shows a transition (or crossover) between a relaxed ($\bar{\zeta}=0$) regime, where supercoiling exists as twist which annihilates, and a writhed ($\bar{\zeta}>0$) regime at forces $F_{p}^{*}\simeq 20\epsilon/\sigma$ (Fig.4), where plectonemes form.

In line with Ref. [13], we can estimate the extent of residual positive or negative supercoiling induced by an elongating RNAP as $vL/D_{\rm Tw}$, where $v$ is the polymerase velocity, and we use the twist diffusion coefficient as initially there is no writhe. Following Ref. [35], given the torsional rigidity of DNA, plectonemes should form when the supercoiling density exceeds $\simeq 0.01$. This criterion suggests plectoneme formation should start from $v\sim 0.005$ bp/$\tau_{Br}$, or $F_{p}^{*}\simeq 50\,\epsilon/\sigma$, roughly in agreement with the value found numerically. It also predicts that plectoneme formation should require a smaller $v$ (hence smaller $F_{p}$) for larger $L$, as we find numerically (Fig. 4). Whilst our simple argument works reasonably well, it disregards the conformational dynamics of the polymer, and in particular the diffusion of tension along the molecule. Estimating the latter as $\sim\sigma^{2}/\tau_{Br}\simeq 10$ bp2/$\tau_{\mathrm{Br}}$ we expect that this contribution may quantitatively affect our estimate for residual supercoiling; it will also introduce an asymmetry between plectonemes formed ahead and behind of the RNAP (as in Fig. 1).

Discussion

In this work we used a single nucleotide resolution coarse-grained model for dsDNA to study the non-equilibrium generation and dynamics of supercoiling by RNA polymerase. Importantly, because our simulations can resolve the full 3-D supercoiling dynamics, they go beyond recent works that have studied this phenomena in 1-D [28, 14]. Our simulations confirm that when the rotational motion of the polymerase is hindered, its transcriptional activity generates TSDs, where DNA in front of the RNAP becomes positively supercoiled, and the DNA behind it becomes negatively supercoiled [11].

The first novel result of our work is that RNAP elongation can effect an “action-at-a-distance”. Contrary to typical textbook pictures [11], we found that the transcriptionally generated twist quickly diffuses away from the RNAP and generates writhe (in the form of plectonemes) far away from it. When the force exerted by the RNAP is not enough to break the H-bonds in the DNA double helix, the RNAP moves in a stick-slip fashion and the position of the plectonemes is halfway to the closest topological barrier. When the force is larger, there is symmetry breaking due to DNA tension, which disfavours plectoneme formation ahead of the RNAP (and favours it behind it, Fig. 1). This “action-at-a-distance” can destabilise nucleosomes, leading to DNA unwrapping long before there is any physical contact between the RNAP and the proteins (Fig. 2). This result is likely relevant in vivo for eukaryotes [30, 36], since it implies that transcription itself can de-stabilise downstream nucleosomes, thereby facilitating further RNAP progress.

The second key finding is that there are two dramatically different modes of supercoiling relaxation. We observed that twist rapidly diffuses along DNA and its diffusion constant is about $D_{\rm Tw}\simeq 500$ bp2/$\tau_{\mathrm{Br}}$. To the contrary, writhe, in the form of plectonemes, relaxes much more slowly and its effective diffusion coefficient is only about $D_{\rm Wr}\simeq 1$ bp2/$\tau_{\mathrm{Br}}$ (consistent with writhe relaxation requiring large-scale rearrangements of the DNA). Additionally, writhe relaxation entails at least two characteristic time scales which are associated with the dissipation of internal stress localised at plectonemes, followed by slower relaxation of delocalised writhe. The presence of two relaxation times is reminiscent of the behaviour of colloidal glasses [37], and it would be intriguing to analyse this parallel further in the future. Finally, for both twist and writhe, positive and negative supercoils relax differently, consistent with the asymmetry in deforming a chiral molecule such as DNA beyond the linear regime [35, 27].

Our results may be tested with in vitro experiments with multiple polymerases on plasmids of different size, or with tethered linear DNA in “curtain” arrangements [38]. Topological barriers can be introduced by including DNA-binding proteins which restrict DNA twist, or non-elongating RNAPs. It would also be useful to analyse the relaxation of plectonemic supercoiling and the anomalous diffusion of writhe via imaging techniques such as those used in Ref. [38]. By performing experiments on reconstituted chromatin fibres one could finally examine whether and when nucleosomes are destabilised or evicted downstream of the advancing RNAP [30].

In the future, we suggest it would be desirable to extend this modelling to include the activity of topological enzymes such as topoisomerases [2] and to account for the asymmetry between the major and minor groove, giving rise to a direct coupling between twisting and bending [39, 40] which is not included in our model. Finally a new direction to take our work further would be the study of transcriptionally-induced supercoiling on plasmids with multiple genes and topological barriers. This set-up could yield a 3-D gene network with purely mechanical regulatory interactions.

Materials and Methods

We simulate dsDNA loops using the model described in Ref. [16]. Each nucleotide is represented by a rigid body formed by a bead and a patch. Beads are connected into a chain, and two chains wind round each other to form a double helix. We use the LAMMPS molecular dynamics software [41] to perform Brownian dynamics simulations where the position of each nucleotide is determined by a Langevin equation (including terms for inter-nucleotide interactions and an implicit solvent). Full details of the models and how we calculate the various quantities (twist, writhe, local and unsigned writhe etc.) are given in SI Appendix.

Acknowledgement

This work was supported by ERC (CoG 648050, THREEDCELLPHYSICS). Y. A. G. F. acknowledges support from CONACyT PhD Grant 384582.

Supporting Information Text

1 A Nucleotide Resolution dsDNA Model

We  performed Langevin dynamics simulations using the LAMMPS software [41], and employed the model described in Ref. [16]. In that reference the reader can find a detailed explanation of its parametrization and validation. For completeness, here we summarize the fundamental characteristics of the model.

A single nucleotide is represented by a rigid body made of two particles: a bead with an excluded volume, diameter $\sigma=1$ nm, and a patch located at the edge of the bead and with no excluded volume (see Fig. S1a). The centre of the bead represents the DNA sugar-phosphate backbone and the patch represents the nitrogenous base. The excluded volume is assigned through a truncated and shifted Lennard-Jones (LJ) potential,

[TABLE]

where $r$ is the separation between bead centres and $\epsilon=k_{B}T$ is an energy scale.

We used FENE bonds of equilibrium length $d_{bp}=0.46$ nm to connect adjacent beads to create a chain that resembles a single-stranded (ss) DNA. The potential has the form

[TABLE]

where we set the maximum bond length $R_{0}=0.68~{}\sigma$, and the energy $\epsilon_{\rm bb}=30~{}k_{B}T$.

The stacking interaction between bases is set by a Morse potential between consecutive patches in the same strand, with equilibrium distance $r_{0\text{, stack}}=0.34$ nm. This potential reads

[TABLE]

with energy $\epsilon_{m}=30~{}k_{B}T$, and $\lambda=8$ is a parameter related to the width of the potential.

In addition, the planarity of these bases is imposed through a harmonic potential

[TABLE]

which constrains the angle $\alpha$ formed between two consecutive patches and one bead (all in the same strand) to have an equilibrium value $\alpha_{0}=90\degree$. In this case $\epsilon_{p}=200~{}k_{B}T$.

The ratio of the backbone equilibrium distance ($d_{bb}$) and the stacking distance($r_{0\text{, stack}}$) imposes a local twist of approximately $36\degree$. However, this is not enough to set the right-handedness of the molecule. Therefore, in order to model this feature we used a dihedral interaction between the particles of two consecutive nucleotides in the same strand, shown in Fig. S1b and given by

[TABLE]

where $\phi$ is the dihedral angle between planes formed by the quadruplets of monomers in two consecutive nucleotides, and $d=180\degree-36\degree$ is a phase angle related to the equilibrium helical pitch. This potential also allows us to regulate the torsional stiffness of the molecule, by varying the magnitude of $\epsilon_{d}$. As detailed in Ref. [16] we set this to give realistic parameters for DNA.

In Fig. S1c, two complementary nucleotides are held together by a breakable harmonic spring which mimics the behaviour of the hydrogen bond (HB). This is given by

[TABLE]

where $r_{0\text{,hb}}=0$ is the equilibrium bond distance, $r_{c\text{,hb}}=0.3$ nm is the cut-off distance beyond which the bond breaks ($U_{hb}(\rm r)=0$) and $\epsilon_{\rm HB}$ is the strength of the HB – we typically set this to $6$ $k_{B}T$, but also consider different values in selected cases as detailed in the main text and below.

The stiffness of the DNA molecule is controlled by means of a Kratky-Porod potential. This potential regulates the angle $\theta$ between three consecutive bases along the same strand and is given by

[TABLE]

The energy $\epsilon_{b}$ is tuned to reproduce the persistence length of DNA. The values used for all parameters and the potentials are fully described in Ref.[16].

If we try to assemble the nucleotides described above into a dsDNA molecule in its B form, the spherical shape of the nucleotides would produce a large steric repulsion between base-pairs (bp) - i.e., the geometry of DNA is not compatible with spherical bodies of this size as building blocks. The solution to this issue in our model was to use two types of beads along the same strand: sterically interacting beads (shown as small solid blue spheres for one strand and red for the other in Fig. S1d) are intercalated by two “ghost” beads (represented as small grey spheres). While theses grey beads do not interact sterically along the same strand, they do interact with all beads on the complementary strand with an excluded volume of $\sigma/2$. This choice ensures that the two strands will not cross each other.

Finally, the position $\mathbf{r}$ of the nucleotides in the system obeys the Langevin-Equation

[TABLE]

where $\mathcal{U}$ is the total potential field experienced by the nucleotide, $m$ is the mass of the nucleotide, $\xi$ is the friction and $\mathbf{\eta}$ is a white noise term with zero mean and satisfying the fluctuation dissipation theorem

[TABLE]

along each Cartesian component (denoted by Greek letters). A similar Langevin equation describes the orientation of each patch-bead rigid body. In LAMMPS, integration of these equations is performed using a velocity-Verlet algorithm [41].

Above we have defined the simulation length and time units $\sigma=1$ nm and $\epsilon=k_{B}T$ respectively, and we use these throughout this work. The unit simulation unit for force can then be expressed as $\epsilon/\sigma$. A natural simulation time unit is the time it would take a bead of diameter $\sigma$ to diffuse across its own diameter; this is known as the Brownian time $\tau_{\rm Br}$ and is related to the diffusion constant for a bead through $D=\sigma^{2}/\tau_{\rm Br}$ (which itself is related to the friction though the Einstein relation $D=k_{B}T/\xi$). To map these to physical units one must set a temperature and solvent viscosity (which can be related to $\xi$ e.g. through Stokes’ law); since this will differ for e.g. an in vitro and an in vivo system we keep our results general and use the simulation units.

2 RNA Polymerase Model

The RNA polymerase (RNAP) is modelled as a rigid body with two parts (see Fig. S2(a)): a ring (shown in black) and a segment across its diameter (shown in yellow) which is referred to as the “crossbar” in the main text. The crossbar is placed in between two consecutive bps in the initial configuration (Fig. S2(b)) with its axis perpendicular to the DNA centre line. It is responsible for breaking the hydrogen bonds and it is made up from seven small beads. The six beads coloured in yellow have an excluded volume of 0.34 nm to repel and break the HB of the dsDNA; the bead in the centre of the crossbar (coloured in green) can also exert a force ($\mathbf{F}_{p}$) on the nearest nucleotides (those within a distance of $0.5$ nm) which is directed perpendicular to the crossbar. Since in our simulations we have fixed the polymerase in place in 3-D space, the applied force leads to motion of the DNA through the ring (in a clockwise direction when $\mathbf{F}_{p}>0$ and in a counter-clockwise direction otherwise). Finally, the ring is made by ten beads ($1$ nm in diameter each) surrounding the crossbar and preventing the RNAP complex from being expelled from the DNA molecule at any point during the simulation (Fig. S2(c)). Whilst the RNAP model we use is topologically realistic, we note that its geometric size is smaller than in reality as, for instance, a bacterial RNAP can be modelled as a $9$ nm $\times$ $9.5$ nm $\times$ $16$ nm ellipsoid [42].

3 Moving Polymerase

In the main text we describe simulations where the position and orientation of the RNAP complex is fixed. In this section we present results from simulations with a mobile RNAP. In this case, the generation and dynamics of supercoiling will also depend on the relative rotational motion between the DNA and RNAP [14]. If, for example, the polymerase crossbar rotates 360° (with the same handedness of the DNA) as it covers a distance of 10 bp along the DNA axis, no supercoiling will be generated, even in the presence of a topological barrier. Therefore, twin supercoiling domain formation requires the rotation rate of the polymerase to be small, meaning that the rotational drag it experiences must be sufficiently large. In vivo this is likely to be the case: a transcribing RNAP is associated with a large number of co-factors [21, 20], there is a nascent RNA connected to the polymerase, and the whole complex may be connected to distant DNA regions or other cellular components.

In our simulations each rigid body (each DNA bead-patch pair or the RNAP) experiences both a translational drag $\xi$ (see Eq. (8) above) and a rotational drag $\xi_{r}$. This drag can be related to an inertial time scale $\tau_{\rm in}=m/\xi$, the time over which the velocity decorrelates (directional information is lost). Here $m$ is the total mass of the rigid body. Similarly for the rotational drag $\tau_{\rm in~{}rot}=I/\xi_{r}$, where $I$ is the moment of inertia. For a sphere, Stokes’ law relates the drag to the viscosity of the solvent $\nu$: $\xi=3\pi\nu\sigma$ and $\xi_{r}=\pi\nu\sigma^{3}$ where $\sigma$ is the diameter of the sphere. The moment of inertia for a solid sphere is $I=m\sigma^{2}/10$. So if we take the simple approximation that the polymerase rotates like a solid sphere, then in simulations where the polymerase can move we should set

[TABLE]

In this case, the rotational drag is small enough such that the RNAP rotates as it moves and no supercoiling is generated. If we instead increase the rotational drag to account for the possible effects listed above, then the RNAP does indeed generate a twin supercoiling domain (Fig. S3 - here we set $\xi_{r}$ so that the RNAP rotates like a solid sphere of diameter $30$ nm).

4 Polymerase-DNA Relative Velocity

Returning to the case of an RNAP fixed in space, in Fig. S6a we show the bp position of the RNAP on the DNA, $x_{p}$, as a function of time for different forces, for simulations with the set-up as in Fig. 1 in the main text. The slope of these curves represents the velocity of the polymerase ($v$), in units of bp$/\tau_{\mathrm{Br}}$. Figure S6b shows this velocity (averaged over $10$ simulations) as a function of $F_{p}$.

The plots show that the RNAP requires a critical force $F^{*}$ to break up the H bonds and move steadily along the DNA backbone. For $F_{p}<F^{*}$, motion is slow and depends on thermally activated breakage of the H bonds. For forces sufficiently larger than $F^{*}$ the velocity of the RNAP grows linearly with $F_{p}$. The crossover between slow and fast motion is similar to that exhibited by the Prandtl-Tomlinson model for frictional motion at the nano-scale [25]. In line with this interpretation, we find that the critical force $F^{*}$ separating slow and fast RNAP motion increases with the H bond energy $\epsilon_{\rm HB}$ between complementary nucleotides. In Fig. S6b we show the results for three different H bond energies $3,6,$ and $12$ $k_{B}T$ (set in such a way to prevent the DNA melting in the negatively supercoiled region). These energies give rise to $F^{*}=20,25,$ and $29~{}\epsilon/\sigma$ respectively.

5 Computing the Local and Total Twist

To obtain the local twist angle $\theta(n)$ between base pairs $n$ and $n+1$ we followed the procedure described in Ref. [43]. We consider the circular DNA molecule as a discrete chain consisting of a set of $N$ vertices, with their locations given by $\mathbf{r}(n)$ for $n=1,\cdots,N$. In our model, $\mathbf{r}(n)$ corresponds to the mid-point of the line connecting the centres of the two beads which form the $n$-th bp. We then define a local reference frame at each vertex. This is specified by the unit tangent vector $\mathbf{\hat{t}}(n)=(\mathbf{r}(n+1)-\mathbf{r}(n))/\lvert\mathbf{r}(n+1)-\mathbf{r}(n)\rvert$, a unit normal vector $\mathbf{\hat{f}}(n)$ obtained from the projection of the vector connecting the two beads in a bp onto the plane perpendicular to $\mathbf{\hat{t}}(n)$, and a third vector $\mathbf{\hat{m}}(n)=\mathbf{\hat{t}}(n)\times\mathbf{\hat{f}}(n)$. Note that since the system forms a ring, the last tangent vector, $\mathbf{\hat{t}}(N)$, joins the last bp to the first.

We define the binormal vector $\mathbf{\hat{b}}(n)=\mathbf{\hat{t}}(n-1)\times\mathbf{\hat{t}}(n)/\lvert\mathbf{t}(n-1)\times\mathbf{t}(n)\rvert$ and the turning angle $\alpha(n)=\arccos(\mathbf{t}(n-1)\cdot\mathbf{t}(n))$ between the incoming and outgoing tangent vectors. For the case $n=1$, the vector $\mathbf{t}(n-1)$ refers to the last tangent vector.

Within this framework, we can rotate the vectors $\mathbf{m}(n-1)$ around $\mathbf{b}(n)$ by an angle $\alpha(n)$, to find the material-shifted vector $\mathbf{m}^{*}(n)$. Finally, the local twist angle can be computed with the use of the following relation,

[TABLE]

The total twist $\mathrm{Tw}$ is found from adding the value of the local twist along all the bps

[TABLE]

6 Computing Writhe, Unsigned-Writhe and Linking Number

The writhe of a closed curve $C$ is defined through the Gauss double integral,

[TABLE]

where $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are any two points along the curve. In the discrete case, the curve $C$ is represented by the set of vectors joining consecutive bps ($\mathbf{t}_{n}$, defined as in the previous section) and the double integral becomes a double sum. The numerical evaluation of this integral can be computed following the methods presented in Ref. [33], where the authors show that Eq. (13) along a curve with $N$ segments can be discretised as

[TABLE]

where $\Omega_{ij}/4\pi$ is the Gaussian integral along the segment $i,j$. These are denoted by $\mathbf{t}_{i,i+1}$ and $\mathbf{t}_{j,j+1}$, which represent two different vectors linking any two consecutive vertices ($i$-$i+1$ and $j$-$j+1$) along the DNA axis. In this case, the absolute value of the Gauss integral multiplied by $4\pi$ represents all the perspectives in which the two previous vectors appear to cross in the solid angle $\Omega^{*}_{i,j}$. This quantity can be found through the definition of the unitary vectors

[TABLE]

and then

[TABLE]

The contribution to the Gauss integral depends on the crossing orientation between two segments and therefore

[TABLE]

The total writhe is obtained by plugging this result into Eq. (14). The unsigned writhe ($\bar{\zeta}$), which we refer to in the main text, is computed by neglecting the sign function in the previous equation before adding all the contributions of the segments $i$ and $j$. Finally, once the values of both writhe and twist have been found, the linking number is computed by adding these quantities, $\mathrm{Lk}=\mathrm{Tw}+\mathrm{Wr}$.

An example of the results of the calculations of twist global twist, writhe and linking number versus time is shown in Fig. S4, for the dynamics corresponding to the simulation procedure outlined in Fig. 1 in the main text. It can be seen that the linking number is conserved at all times, as is topologically required for a ring-like DNA molecule.

7 Length of the First Plectoneme and Distance to the Polymerase

In order to calculate the length ($l$) of a plectoneme, and its distance ($d$) from the RNAP, we examine the contacts made between non-adjacent bp – where a contact is defined when the separation between two bp is less than a threshold distance $r_{\rm threshold}$. A contact map is a visual representation of all the contacts within the dsDNA molecule, where a matrix is obtained by filling the $m$-$n$ element with 1 if two bp are in contact, and 0 otherwise. Self-contacts are not allowed and therefore values on the matrix diagonal are zero. We set $r_{\rm threshold}=3\sigma=3$ nm; since this value is larger than the rise of the DNA, we can expect that close to the diagonal of the matrix the elements will have a value of one.

As in Sec. 5, the middle point between two complementary beads represents the position $\mathbf{r}_{n}$ of the $n$-th bp in our model. In Fig. S5a we show schematically this position for a small section of the DNA molecule, lying between bp $i$ and bp $j$. If this piece of DNA forms a plectoneme like the one shown in Fig. S5b, that would mean that the contour length ($l$) of the plectoneme is $\lvert i-j\rvert$ (in bps). With a threshold choice of $r_{\rm threshold}=3\sigma$, as well as beads $i$ and $j$, beads $i+1$ and $i-1$ and beads $i+2$ and $j-2$, and so on to bead $(i+j)/2$, also tend to be in contact. In the contact map, the signature of this configuration is a line of ones perpendicular to the diagonal of the matrix.

In Figs. S5c-d, we show the contact map after the index of the bps have been shifted so the RNAP is always located in the middle of the molecule (at position $L/2=500$). Panel (c) represents a typical case of the first plectoneme appearing in our simulations with negative writhe (behind the polymerase), while (d) represents a typical case with positive writhe (in front of the polymerase). The distance to between the RNAP and a plectoneme is then either $\lvert N/2-i\rvert$ or $\lvert N/2-j\lvert$, whichever is smaller. We generate a contact map every $\delta t$ time steps, and take the measurements for $l=\lvert i-j\rvert$ and $d=\min\left[\lvert N/2-i\rvert,\lvert N/2-j\lvert\right]$ at the first time point where a plectoneme signature appears.

8 Dependence of Nucleosomal Wrapping on Supercoiling

As described in the main text, we model a histone octamer as a large bead ($10$ nm in diameter) with $21$ DNA-binding patches on it, following a left-handed helical path. The patches have a short range attractive interaction with the DNA beads. In addition, the separation between these patches is set so that overall they favour wrapping of $150$ bp of DNA ($1.7$ turns around the histone octamer), as in real nucleosomes. We observed in our simulations that a $1000$ bp long dsDNA molecule, initialized with $Lk=Lk_{0}$ (relaxed state), can fully wrap around the histone octamer, without causing plectonemic writhe formation at any other location (Fig. S7a, with DNA bp indices shown in black, and the index of the histone octamer patches in red).

An initial “set-up” simulation is run starting from a configuration where the nucleosome bead is placed adjacent to a relaxed ($Lk=Lk_{0}$) DNA molecule with an inactive RNAP; during this simulation the DNA wraps around the bead and the nucleosome self assembles. Subsequent simulations use this wrapped configuration as a starting point; a TB is positioned bps $230$–$240$, and the RNAP is activated, as detailed in the main text.

In Fig. 2 of the main text we show bar plots indicating mean values of the local twist at different regions, as well as kymographs showing the (un)wrapping behaviour of the nucleosome when subject to supercoiling of different signs. In Fig. S7 we plot the full local twist profiles for every bp in the DNA loop. Figure S7a shows the profile at the beginning of the simulation when the nucleosome is fully formed, while Figs. S7b and c show the profile after the RNAP has moved for the cases where the nucleosome is subjected to positive and negative supercoiling respectively. Each plot is an average over 50 independent simulations.

9 Twist Diffusion Coefficient

If the contour length of the DNA ring is expressed through the continuum variable $x\in$ $[0,L]$, with $L$ the number of bps along the molecule, the diffusion equation of the local twist $\theta(x,t)$ at time $t$ and position $x$ is then given by

[TABLE]

where $D_{\rm Tw}$ is the twist diffusion coefficient. We then use an initial condition where a region located between points $x_{i}$ and $x_{f}$ ($x_{i}<x_{f}$) has a local twist $\theta$, different from that of B-DNA in its relaxed form ($\theta_{0}=36\degree$). Then, $\Delta\theta=\theta-\theta_{0}$ represents the excess (or deficit) of twist in this region, and the initial condition can be expressed as:

[TABLE]

where $H$ is the Heaviside step function. The solution of Eq. (21) with the initial condition Eq. (22) is

[TABLE]

where $\text{erf}(x)$ is the error function.

Finally, the data obtained from a simulation at any time step $t$, can be fitted to Eq. (23) to obtain the twist diffusion coefficient.

10 Local Writhe Field

In order to study the relaxation of plectonemes in our model, we ran simulations starting from a configuration of a DNA molecule where twin-supercoiling domains have been generated by the action of an RNAP (as in the right-hand snapshot of Fig. 1b in the main text for example). We then remove both the topological barrier and the RNAP complex, and we analyse the evolution of this new system to its relaxed state. Under these conditions, the two opposite sign plectonemes will relax and eventually disappear. One way of studying this relaxation process is by defining an effective local writhe field $\chi(n,t)$ at position $n$ along the DNA and time $t$, that can be tracked directly from our Brownian dynamics (BD) simulations. To define $\chi(n,t)$ we consider a segment of the DNA as a “sliding window” centred at bp $n$. The length of this window should be larger than or equal to the length of the largest plectoneme in our initial configuration. We then create an artificially closed curve that first follows the open path of the DNA segment contained within the window, and then is extended by adding extra beads along the tangent direction at both ends of this path. After reaching a sufficient distance from the plectoneme (to minimise the chance of introducing writhe artificially via the closure procedure), the curve can be closed by defining a vector joining the two new terminal beads. We can then we compute the writhe contained in the resulting ring by solving the Gaussian integral (Eq.(13)). Next, we move the centre of the window by one bp and repeat the procedure; this continues until we have covered the whole contour length of the DNA. In this way we obtain $\chi$ as a function of $n$ for a given time point. The temporal evolution of the field $\chi$ is followed by repeating the calculation at different time-steps. We repeated the procedure mentioned above for four cases (case 1 is shown in Fig. 3b in the main text, and cases 2-4 correspond to Figs.S8a-c), which differ in the number of plectonemes in the initial configuration and the magnitude of $\chi$ stored in these plectonemes.

The writhe field plots shown in Figs. S8a-c reveal how the amplitude of the local writhe in the plectonemes decreases as the system evolves towards the relaxed state. In the initial configuration, the maximum of $\chi(n,0)$ corresponds to the writhe stored in the plectoneme located in front of the RNAP, whereas the minimum corresponds to the writhe stored in the plectoneme behind the RNAP complex. The position of the two extremes in this plot coincide with the indices of the bps at the tip of the plectonemes.

11 Writhe Relaxation Dynamics

In Fig. S9 we show the temporal evolution of the maximum of $\chi(n,t)$ for the cases mentioned in the previous section. We observe that in all cases the relaxation of the local writhe involves two timescales: an initial rapid relaxation where the magnitude of the writhe decreases as $e^{-t/\tau_{1}^{+}}$, and then the process slows down, decreasing as $e^{-t/\tau_{2}^{+}}$ (with $\tau_{2}^{+}>\tau_{1}^{+}>0$). Similarly, the evolution of the magnitude of the minimum of $\chi(n,t)$ also shows two timescales $\tau_{1}^{-}>\tau_{2}^{-}$ (Fig. S9 inset).

Values of the decay constants measured for each case simulated are given in Table S1. We find that in each case $\tau_{1}^{+}<\tau_{1}^{-}$ (except in case 1, where the plectoneme with positive writhe is so short that the first time-scale cannot be computed) and $\tau_{2}^{+}<\tau_{2}^{-}$. This implies that positive writhe relaxes more quickly than negative writhe, in line with our observation that the formation of the first negative plectoneme occurs earlier than the first positive plectoneme in our simulations (see main text). Put another way, if negative twist is energetically more costly than positive, this explains why the negative side writhes first in the simulation where supercoils are generated; in the simulation where supercoils relax, the fact that negative writhing relaxes more slowly implies that the relaxation pathway involves conversion to (more costly) negative twist.

12 Bending Energy During Writhe Relaxation

In the previous section we showed that the relaxation of writhe involves at least two time-scales. One explanation for this behaviour is that the tension accumulated in the relative short plectonemes drives the fast relaxation at the beginning. When a large part of this tension has been released by decreasing the local writhe in the plectoneme, the relaxation slows down accounting for the second slower timescale. In order to corroborate this hypothesis, we now analyse the evolution of the bending energy ($E_{\mathrm{bend}}$) of the DNA molecule during the overall relaxation process.

In our model the energy associated with bending of the dsDNA can be calculated as

[TABLE]

where $\mathbf{t}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i+1}/|\mathbf{r}_{i}-\mathbf{r}_{i+1}|$ is a tangent vector at bp $i$, and $l_{p}=50$ nm and $l=0.34N$ nm are the DNA persistence length and bp length respectively. At different times the shape of the molecule will change and the temporal evolution of the bending energy can be tracked. After averaging over five different trajectories we obtain the curves shown in Fig. S10, each one representing the evolution from different initial configurations (cases 1-4 in the previous section). Clearly, the configurations with larger plectonemes at $t=0$, start with the highest bending energy. We see that the value of $E_{\mathrm{bend}}$ decays quickly at the beginning of the relaxation and then stabilizes, presumably corresponding to the two relaxation timescales observed in the writhe dynamics.

13 A Model for Writhe Relaxation

Here we discuss a model for the relaxation of writhe. As seen in the main text, while the twist appears to diffuse at a constant rate, the writhe displays two distinct relaxation regimes with a first fast decay and a subsequent slower dynamics.

We describe the decay of the maximum $w(t)=\max[\chi(x,t)]$ (or $v(t)=\min[\chi(x,t)]$) as reported in the main text (Fig. 3c and Fig. S9) using a phenomenological ODE,

[TABLE]

with

[TABLE]

The chosen functional form connects the early time diffusive relaxation (characterised by a relaxation constant $d_{1}$) with the late time regime (characterised by a constant $d_{2}$) through a smooth sigmoidal function parametrised by a thickness $\sim b^{-1}$ and a transition point $w_{c}$.

We numerically solve Eq. (25) for different choices of the four free parameters $d_{1},d_{2},b$ and $w_{c}$ and minimise the sum of the square residuals $R^{2}=\int_{0}^{T}dt\left[w(t)-w_{o}(t)\right]^{2}$, where $w_{o}$ is the maximum of the writhe computed from our BD simulations. The result from this procedure is plotted in Fig. S11 where we show $\tilde{\chi}(t)=w(t)/w(0)$ (and $v(t)/v(0)$) as a function of time. As pointed out in the main text, all curves display similar behaviours despite starting from different values of local writhe. The transition in between the two regimes is marked by the critical writhe $w_{c}$ which is found to roughly be a fixed fraction of the initial value of writhe – i.e. $w_{c}/w(0)\simeq 0.66$.

Case Maximum[$\chi(n,t)$] -Minimum[$\chi(n,t)$]

$\tau_{1}^{+}$ [$\tau_{\mathrm{Br}}$] $\tau_{2}^{+}$ [$\tau_{\mathrm{Br}}$] $\tau_{1}^{-}$ [$\tau_{\mathrm{Br}}$] $\tau_{2}^{-}$ [$\tau_{\mathrm{Br}}$]

1 $---$ 27213$\pm$623 4915$\pm$207 26036$\pm$416

2 2580$\pm$120 20443$\pm$373 3469$\pm$79 22811$\pm$405

3 3066$\pm$119 27976$\pm$284 4531$\pm$206 68989$\pm$959

4 2991$\pm$131 38252$\pm$333 4857$\pm$161 59976$\pm$488

average 2879$\pm$71 28471$\pm$212 4443$\pm$86 44453$\pm$306