Plectoneme dynamics and statistics in braided polymers

TL;DR

This paper investigates the formation, dynamics, and statistics of plectonemes in braided polymers using simulations and analytics, revealing phase transitions, kinetic behaviors, and a stochastic theory for twist-writhe coupling.

Contribution

It introduces a new order parameter and a stochastic field theory to explain plectoneme behavior and phase transitions in braided polymers.

Findings

Identification of phase boundary via plectonemic fraction

Observation of phase separation kinetics and plectonemic hopping

Development of a stochastic model explaining coarsening laws

Abstract

Braids composed of two interwoven polymer chains exhibit a buckling transition whose origin has been explained through the onset of plectonemic structures. Here we study, by a combination of simulation and analytics, the dynamics of plectoneme formation and their statistics in steady state. The introduction of an order parameter, the plectonemic fraction, allows us to map out the phase boundary between the straight braid phase and the plectonemic one. We then monitor the formation and the growth of plectonemes, observing events typical of phase separation kinetics for liquid-gas systems (fusion, fission, 1D Ostwald ripening), but also of DNA supercoiling dynamics (plectonemic hopping). Finally, we propose a stochastic field theory for the coupled dynamics of twist and local writhe which explains the phenomenology found with Brownian dynamics simulations as well the power laws underlying the coarsening of plectonemes.

Figures

Equations

$$U_{\text{mbox}FENE}=-\frac{K_{\text{mbox}f}}{2}\sum _i^{N-1}{\left(\frac{R_0}{\sigma }\right)}^2\ln \left[1-{\left(\frac{r_{i,i+1}}{R_0}\right)}^2\right],$$

$$\frac{U_{\text{mbox}LJ}}{4ϵ}=\sum _{i,j>i}^N\left[{\left(\frac{\sigma }{r_{i,j}}\right)}^{12}-{\left(\frac{\sigma }{r_{i,j}}\right)}^6+\frac{1}{4}\right]\theta (2^{\frac{1}{6}}\sigma -r_{i,j}),$$

$$U_{\text{mbox}BEND}=\sum _{i=2}^{N-1}{K}_{\text{mbox}b}\left[1+\cos {\varphi }_i\right],$$

$$\frac{U_{\text{mbox}LJ}^{\prime }}{4ϵ}=\sum _{i,j}^N\left[{\left(\frac{\sigma }{r_{i,j}}\right)}^{12}-{\left(\frac{\sigma }{r_{i,j}}\right)}^6+\frac{1}{4}\right]\theta (2^{\frac{1}{6}}\sigma -r_{i,j}),$$

$$U_{\text{mbox}bot}=K_{\text{mbox}w}(r-r_c{)}^2,$$

$$\leavevmode\resizebox{422.77661pt}{}{$\frac{U_{\mbox{\scriptsize{top}}}}{4\epsilon}=\sum_{\{b,w\}}\left[\left(\frac{\sigma}{d_{b,w}-\Delta}\right)^{12}-\left(\frac{\sigma}{d_{b,w}-\Delta}\right)^{6}+\frac{1}{4}\right]\theta(2^{\frac{1}{6}}(\sigma+\Delta)-d_{b,w})\;,$}.$$

$$m{\ddot{\mathrm{r}}}_i=-\gamma {\dot{\mathrm{r}}}_i+{\xi }_i+{\nabla }_{{\mathrm{r}}_i}U({\mathrm{r}}_i),$$

$$⟨{\xi }_{i,\alpha }(t)⟩$$

$$⟨{\xi }_{i,\alpha }(t)\cdot {\xi }_{j,\beta }(t^{\prime })⟩$$

$$\frac{\partial \mathrm{Tw}}{\partial t}$$

$$\frac{\partial \mathrm{Wr}}{\partial t}$$

$$g(\mathrm{Tw},\mathrm{Wr})$$

$$\frac{\partial \mathrm{Tw}}{\partial t}$$

$$\frac{\partial \mathrm{Wr}}{\partial t}$$

$$g(\mathrm{Tw},\mathrm{Wr})$$

$$⟨\eta (x,t)\eta (x^{\prime },t^{\prime })⟩=\delta (x-x^{\prime })\delta (t-t^{\prime }),$$

$$\frac{\partial \left(\mathrm{Tw}+\mathrm{Wr}\right)}{\partial t}=\frac{\partial }{\partial x}\left[D_t\frac{\partial \mathrm{Tw}}{\partial x}+D_w\frac{\partial \mathrm{Wr}}{\partial x}\right].$$

$$\frac{\partial \sigma }{\partial t}$$

$$\frac{\partial \delta }{\partial t}$$

$$\delta =\frac{\beta -\alpha -\gamma \Theta \left(|\sigma |-{\sigma }_0\right)}{\beta +\alpha -\gamma \Theta \left(|\sigma |-{\sigma }_0\right)}\sigma .$$

$$\frac{\partial \sigma }{\partial t}=M\frac{{\partial }^2{\mu }_{\mathrm{eff}}}{\partial x^2}\equiv M\frac{{\partial }^2}{\partial x^2}\left(\frac{\partial f_{\mathrm{eff}}}{\partial \sigma }\right),$$

$$M{f}_{\mathrm{eff}}$$

Full text

I Introduction

Polymer braids lie at a crossroad between topology, statistical mechanics, and biological physics. Within cells, two DNA double-helices inevitably intertwine to form a braided molecule, or catenane, when replication terminates Mariezcurrena and Uhlmann (2017). Braids also arise in biology with interwoven collagen helices Orgel et al. (2006), or amyloid fibrils Ionescu-Zanetti et al. (1999). The biophysical properties of braided polymers can be quantitatively studied with single molecule experiments in vitro Terekhova et al. (2014); Cejka et al. (2012); Brahmachari et al. (2017). In nanotechnology, semiflexible polymers tethered to colloidal particles may be used to create biomimetic braids at near-molecular scale Goodrich and Brenner (2017).

From a statistical mechanics viewpoint, these systems are of fundamental interest in view of the "buckling" transition between a straight braid phase, where the braid centreline fluctuates around a straight line, and a plectonemic phase, where the centreline writhes in 3D. The plectonemic phase arises for sufficiently large "catenation number" – essentially the linking number between the two polymers in the braid. The thermodynamics of the transition is understood as a competition between torsional stresses stored in a twisted braid and bending energy cost associated with plectoneme formation, which can be formulated at the mean field theory level Brahmachari and Marko (2017); Marko (1997); Neukirch and Marko (2011). The underlying physics is similar to that of buckling in supercoiled DNA under tension, although in braids the absence of H-bonding between the intertwining molecules leads to important qualitative distinctions Brahmachari and Marko (2017).

A quantitative understanding of the transition between straight and plectonemic phase in braided polymers, going beyond mean field, is still elusive, although Monte-Carlo simulations have mapped out the phase diagram for DNA braids under tension, by monitoring the location of the singularity in the plots of braid extension versus catenation number Charvin et al. (2005). The recent experiments in Ref. Brahmachari et al. (2017) also showed evidence of a multimodal distribution of braid extension (end-to-end distance), which arises because the population of braids in equilibrium has multiple plectonemes at large catenation.

In sharp contrast to this body of work on the thermodynamics of the transition between the straight and plectonemic phase, the dynamics in the latter phase has received much less attention. Experimentally, the main study in this field has only addressed the diffusional dynamics of DNA supercoils van Loenhout et al. (2012), whereas theory or simulation studies on braid dynamics are to date altogether lacking. To fill this gap, here we present Brownian dynamics (BD) simulations of two braided semiflexible polymers at fixed values of the catenation number and under a stretching force.

After introducing the model of braided chains under torsion and the relevant observables (see Section II), in section III we first map the boundary between the straight and the plectonemic phase and then we focus on the dynamics of the braid following a quench of variable depth into the latter. Our key finding is that the kinetics of plectonemic growth resembles phase separation in liquid-gas systems. Thus, first twist is converted into writhe to nucleate formation of a plectoneme. Later on, plectonemes coarsen via a combination of fusion and 1D Ostwald ripening. In steady state, we also observe plectonemic “hopping”, whereby a plectoneme seemingly unravels at one place whilst reforming at another, as in the experiments of van Loenhout et al. (2012). Finally in Section IV we propose a stochastic field theory for twist and writhe dynamics which explains the fundamental origin of the phase separation phenomenology we observe as due to the underlying topological conservation of the braid catenation number. Section V is devoted to discussion and conclusions.

II Model and simulation details

Braids are modelled by a pair of semiflexible intertwined chains, each of which is made up by $N$ beads of diameter $\sigma$ and has persistence length $20.6\sigma$ (appropriate, e.g., for double-stranded DNA when $\sigma=2.5$nm) (see Fig. 1). The total energy $U$ of each chain comprises three terms: (i) a finitely extensible nonlinear elastic (FENE) spring potential which accounts for the chain connectivity, (ii) a truncated and shifted Lennard-Jones potential which accounts for excluded volume interactions Kremer and Grest (1990) and (iii) a bending potential which provides the chain with an intrinsic bending stiffness:

[TABLE]

where $r_{i,j}=\left|\boldsymbol{\mathrm{r}}_{i,j}\right|=\left|\boldsymbol{\mathrm{r}}_{j}-\boldsymbol{\mathrm{r}}_{i}\right|$ is the distance between the $i$-th and the $j$-th bead, $R_{0}=1.5\,\sigma$ is the maximum bond length, $\epsilon=k_{B}T$ is the thermal energy of the system, $K_{\mbox{\scriptsize{f}}}=30\frac{\epsilon}{\sigma^{2}}$ is the bond strength, $\theta(x)$ is the Heaviside function and $\phi_{i}$ is the angle between the three consecutive beads $i-1$, $i$ and $i+1$. The value of $K_{\mbox{\scriptsize{bend}}}$ is related, in units of $\sigma$, to the persistence length of the chain $l_{p}\simeq K_{b}/\epsilon$. By using $K_{b}\approx 20\epsilon$ we set $l_{p}=20.6\,\sigma$. Note that beads belonging to different polymers also interact via excluded volume interactions:

[TABLE]

where now $i$ and $j$ are two indices which are running on the first and second chain respectively. To mimic a typical experimental magnetic tweezers set up, each chain is anchored to an impenetrable wall at the bottom while the other pair of ends are fixed to the surface of a top plate that can rotate around the $z$-axis under the action of a constant torque $\tau$. Interactions between the bottom wall and the chains are accounted by the following harmonic potential:

[TABLE]

where $r$ is the distance between the bead and the region surface, $r_{c}=\sigma$ is a cutoff distance below which there is no interaction and $K_{\mbox{\scriptsize{w}}}=200\frac{\epsilon}{\sigma}$ is the spring constant. The top plate is composed by beads (labelled w) of diameter $7\,\sigma$ which interact with the the beads in the two chains (labelled b) through a truncated and shifted Lennard-Jones potential:

[TABLE]

In Eq. 6 the sum runs over all pairs $(b,w)$ with reciprocal distance $d_{b,w}$ and $\Delta=R-\frac{\sigma}{2}$ with $R$ being the radius of the wall bead. Note that in all simulations we set the distance $d$ between the two end beads anchored at the top plate to be $d/L_{0}\simeq 0.17$ where $L_{0}=N\sigma$ is the contour length of each chain Not .

The dynamics of the system is described by the Langevin equation:

[TABLE]

where $\boldsymbol{\mathrm{r}}_{i}$ is the position of the ith bead, $\gamma$ is the friction coefficient and $\boldsymbol{\mathrm{\xi}}_{i}$ is the usual stochastic Brownian noise defined by

[TABLE]

In Eq. 8, $t$ and $t^{\prime}$ are times, $k_{B}$ is the Boltzmann constant, $T$ denotes temperature, $\delta_{i,j}$ is a Kronecker delta, and $\delta(t-t^{\prime})$ a Dirac delta. As characteristic time unit we consider the Brownian time $\tau_{LJ}=\frac{\sigma^{2}}{2D}$ (the time needed by a bead to diffuse across a length equal to its diameter): by using the Einstein relation $D=k_{B}T/\gamma$ and the Stoke’s law $\gamma=3\pi\sigma\eta$ we obtain $\tau_{LJ}=\frac{3\pi}{2}\frac{\sigma^{3}\eta}{k_{B}T}$ where $\eta$ is the viscosity of the fluid inside the flow cell. The time evolution of the system is obtained by integrating numerically Eq. 7 with the LAMMPS software Plimpton (1995) with an integration time step in the range $[0.0025,0.01]\,\tau_{LJ}$.

To simulate the system at equilibrium in the fixed catenation number (or inter-chain linking number) $Ca$ and force $f$ ensemble, initial conditions are prepared by interwining the two chains with the desired value of $Ca$. The latter is monitored by first closing the two chains and then computing their linking number. The system is then relaxed by BD simulations at constant temperature $T$ and reduced force $\tilde{f}=f\frac{\sigma}{k_{B}T}$. Since chains cannot cross the top plate and the bottom wall, the initial catenation $Ca$ is preserved at all times.

A crucial point in our study concerns the detection of plectonemes along the braided filaments (see Fig. 2a). This is performed by first constructing a distance map for each of the two chains. This map is a contour plot of the distances between beads i and j of one chain (e.g. the blue one). [The distance map in Fig. 2a has a cutoff at $r_{cut}=10\,\sigma$  Coronel et al. (2018) as differences in distances of non-contacting monomers are irrelevant.] Plectonemic domains can be identified as darker regions extending perpendicularly from the diagonal in the distance map (Fig. 2b). The portion of these regions closer to and farther from the diagonal represents respectively the apices and the roots of the plectonemes (see Fig. 2(c,d)).

More precisely, the algorithm to detect plectonemes through distance maps works as follows. First, we use a threshold $R_{\mbox{\scriptsize{max}}}$ to identify the set of “contacting beads”, for which $r_{i,j}\leq R_{\mbox{\scriptsize{max}}}$ (these are the darker regions in the maps in Fig. 2b). Then, we measure the positions of the plectoneme extremities (i.e. the points in the ellipses in Fig. 2d) and use these to determine the plectoneme length. In our algorithm, plectonemes retained in the subsequent analysis were required to be larger than an additional threshold $\Delta_{p}$ (measured in units of $\sigma$) Not .

Note that every distance map refers to a single chain: this means that we study plectonemes in each filament separately. Sometimes a chain wraps around the other locally forming a solenoidal structure whose beads satisfy the conditions determining a plectoneme: in this case our method of counting plectonemes can provide $N_{b}$ plectonemes for the blue chain and $N_{r}\neq N_{b}$ plectonemes for the red one. To avoid an overestimate of the number of plectonemes, we assume that this number is given by $\min(N_{b},N_{r})$.

III Results

III.1 Equilibrium phase diagram

The results obtained for different $(Ca,\tilde{f})$ pairs are summarized in Fig. 3. In thermodynamic equilibrium, and over the range of parameters we analyse, the braided polymers can be in one of two phases (Fig. 3c). Either the chains wind around each other without coiling in 3D (straight braid phase, see snapshot in Fig. 3, left), or they end up in a configuration where some of the twist is converted into writhe, leading to a shrinking of the braid along the $z$ direction and the formation of a plectoneme where the backbone of the braid writhes around itself (plectonemic phase, see snapshots in Fig. 3c, right).

To determine the boundary between the straight braid and plectonemic phase, the standard approach is to record the end-to-end distance of the chains along the $z$ direction as a function of $Ca$ – the point at which the slope changes abruptly then corresponds to the transition point (Fig. 3b) Charvin et al. (2005); Brahmachari et al. (2017). Here, we instead look at the statistics of plectonemes by using the algorithm described in Section II and define an overall plectonemic length, $L_{\rm pl}$. The braid is in the plectonemic phase if the order parameter $\langle L_{\rm pl}\rangle/L_{0}$ – the "plectonemic fraction" – is $>0$. The plectonemic fraction as a function of $Ca$ for different values of $\tilde{f}$ is shown in Fig. 3a.

To gain more insight on the nature of the transition we look at the statistics of both the length $\ell_{\rm pl}$ of single plectonemes and their number $N_{\rm pl}$ in steady state. We find that the distribution of plectoneme length is multimodal (Fig. 4a) and that configurations with multiple plectonemes are possible (Fig. 4b). Close to the critical point, plectonemes are short and the length distribution is relatively narrow. Deeper in the plectonemic phase, the distribution widely broadens; additionally plectonemes are longer and more numerous. Multimodality in the length distribution here is much more enhanced and is due to the presence in the population of different classes of configurations, mainly those with a single, two or three plectonemes (Fig. 4b). Interestingly, multimodality in the plectonemic phase was recently reported in a combined experimental and theoretical study, but with respect to distribution of braid end-to-end distance Brahmachari et al. (2017).

III.2 Near-equilibrium dynamics of plectonemes

We now discuss the near-equilibrium dynamics of the system after a quench from the straight braid phase into the plectonemic one by instantly reducing the pulling force. The kinetics depends on the value of $Ca$ and $\tilde{f}$.

Following the quench into the plectonemic phase, the twist soon converts into writhe, with the latter initially localised close to the top plane (where the pulling force is applied). Then plectonemes can grow until they reach their equilibrium size. The kymograph in Fig. 5 monitors the kinetics of plectonemes in the system when $(Ca,\tilde{f})=(31,6)$. In our BD simulations, we observe a number of possible kinetic events. First, plectonemes diffuse slowly along the chain once they are formed. Second, there are fission (Fig. 5 A to B) and fusion (Fig. 5 after B) events, where a plectoneme splits into two, or two close plectonemes merge into one. Third, plectonemes may "hop" from one place to another along the braid, as writhe unravels in one region and reappears in another one (see snapshots from Fig. 5 C to D). Finally, we observe Ostwald ripening events (Fig. 5 E to F) where a small plectoneme is absorbed by a larger one, without touching.

Within the context of DNA supercoiling, plectonemic diffusion was observed in experiments van Loenhout et al. (2012) and simulations Matek et al. (2015), and hopping only in experiments van Loenhout et al. (2012). Inspection of our BD simulations (Suppl. Movie 1) suggests that hopping events are actually preceded by a fluctuation resulting in the moderate shrinkage of an initially large plectoneme. The extra ropelength gained allows nucleation of new small plectonemes, at random positions along the chain, which compete with the original plectoneme for braid length. If one of these outgrow the original plectoneme, the resulting kymograph records a hopping. As hopping requires nucleation, its frequency should strongly depend on noise – this is the case in our simulations and arguably also in experiments van Loenhout et al. (2012). It is then worth noting that fission and fusion events balance in steady state to yield a finite average size for plectoneme. This is similar to liquid-gas phase separation in 1D, where fluctuations inhibit coarsening beyond the correlation length of the system Tailleur and Cates (2008).

The plectoneme dynamics when the quench is deep in the plectonemic phase, compared to the one close to the transition, is very different. This is apparent in Fig. 6 where we contrast the trajectories for quenches to the point $(Ca,\tilde{f})=(32,8)$ (panel (a)) and for $(Ca,\tilde{f})=(32,6)$ (panel (b)). One can notice that both trajectories follow a fast relaxation to equilibration dynamics in which the observable $Z/L_{0}$ (blue line) drops abruptly (this is particularly apparent in panel (b)). This fast decay is due to the sudden reduction of the pulling force at $t=0$, leading to the transition from a straight braid to a braid where plectonemes nucleate.

Close to the plectonemic transition a gas of plectonemes appears: numerous domains nucleate between the lower and the upper fork, but they do not merge as their lifetime is too short. The braid end-to-end distance is weakly influenced by the formation of domains and only occasionally its decrease is clearly visible due to the presence of plectonemes.

Deep in the plectonemic phase instead the system prefers to form longer plectonemes. While the total length of plectonemes is practically constant their number varies. As in Fig. 5, we find events which are typical of liquid-gas phase separation, such as fission (Fig.  6b A to B) and fusion (Fig.  6b C to D). Also in this case the extension of the braid is not strongly influenced by the nucleation of plectonemes and the size of fluctuations is similar to the one seen close to the transition.

III.3 Coarsening dynamics

To characterise more quantitatively the dynamics of plectonemic coarsening following a quench, we analyse in Fig. 7 how the average plectoneme size and number depends on time. The number is non-monotonic with a maximum that is reached at very early times and whose values increase with $Ca$ (see Fig. 7a). This behavior is due to the formation of numerous plectonemes just after the quench and to the following process of coarsening. That the number is larger than $1$ in steady state is consistent with arrested coarsening. The intermediate time behavior is well fitted by a power law $t^{\omega}$, with $\omega\simeq 1/3$ – to a good approximation this is true irrespective of the value of $Ca$ and $\tilde{f}$ (Fig. 7b).

IV Stochastic field theory

The phase separation kinetics between straight and plectonemic phases can rationalized via a simple stochastic field theory, based on two coupled reaction-diffusion equations for the evolution of local twist (${\rm Tw}$) and centreline writhe (${\rm Wr}$), inspired by our BD simulations. These equations read as follows (see appendix A for more details),

[TABLE]

where $D_{t}$ and $D_{w}$ are the diffusion coefficients of twist and (local) writhe respectively, while $\kappa>0$ is a surface tension-like parameter. The reaction terms $\pm g$ are equal and opposite in the two equations to ensure global conservation of the "local catenation" number, ${\sigma}={\rm Tw}+{\rm Wr}$, at all times. The parameters $\alpha$, $\beta$ and $\gamma$ describe interconversion between twist and writhe – we choose $\alpha\geq\beta\geq\gamma>0$ to locally favour formation of a plectoneme for $\sigma>\sigma_{0}$; we also include a noise term, where $\sigma_{\eta}$ denotes noise strength and $\eta(x,t)$ is a Gaussian white noise, uncorrelated in space and time. Eq. 10 (with no noise) is familiar in the literature on Turing pattern formation Murray (2001); Halatek and Frey (2018), where a sufficiently large difference in the diffusivity of the fast and slow fields (here twist and writhe respectively) leads to an instability of the uniform phase (Fig. 8 and Appendix). [In our context, the instability corresponds to phase separation and formation of one or more plectonemes in the braid.] For $D_{t}=D_{w}$, the linking number obeys a simple diffusion equation, hence is uniform in steady state.

Fig. 7c shows the evolution of the average plectoneme size in numerical solutions of our stochastic field theory. Due to the global conservation of $\sigma$ we expect Eq. 10 to be in the same universality class of model B dynamics for liquid-gas phase separation (this is proved in the Appendix), for which we expect a growth exponent $\omega=1/3$ Bray (2002) – as we find numerically in both BD and field theory (Figs. 7b,c). More generally, numerical solutions of our field theory show instances of all dynamical processes seen in simulations. Thus, the kymograph in Fig. 7d shows examples of Ostwald ripening and fusion, whereas in steady state we find hopping and fission if $\sigma_{\eta}$ is sufficiently large (Fig. 9).

V Conclusions

In summary, we have studied the dynamics of polymer braids under tension by numerical simulations. By varying the applied force and the catenation (linking) number between the two chains, we mapped out the transition between the straight braid and plectonemic phase. Our main result is the characterisation of dynamical events and kinetic regimes deep in the plectonemic phase. We find that plectoneme formation shares some features of phase separation in liquid-gas systems and proceeds via a combination of Ostwald-like ripening and fusion events. At late times, there is a balance of merging and breakup of plectonemes which leads to arrest of coarsening with selection of a typical plectoneme size. We also provide evidence of plectonemic hopping, where writhed regions seemingly unravel at one location to appear somewhere else along the braid.

We showed that our numerical data can be explained by a stochastic field theory which models twist and writhe interconversion in our polymer braid. This theory suggests that the growth of plectonemes can be described by a power law, where the exponent is $\sim 1/3$, as befits diffusive phase separation in a liquid-gas system. It also reproduces hopping and concurs with our BD simulations in finding that sufficiently strong noise (i.e., proximity to the transition) is crucial to observe this dynamical mechanism. Both BD and field theory suggest that this phenomenon is due to stochastic nucleation of an additional plectoneme and subsequent competition between the nascent plectoneme and another one in the chain. As the latter unravels, the result is an apparent hopping. In the future, it would be of interest to study quantitatively via single molecule techniques the dynamics and growth laws of plectonemes in DNA and polymer braids, to confirm whether the diffusive liquid-gas phase separation outlined here is confirmed experimentally.

VI Acknowledgments

We thank ERC (CoG 648050 THREEDCELLPHYSICS) for funding.

Appendix

In this Appendix we discuss the details of the continuum reaction-diffusion model whose results have been described in the text. The model consists of two coupled dynamical equations for twist and writhe, which we write down here again for convenience,

[TABLE]

In Eq. 10, $D_{t}$ and $D_{w}$ are the effective diffusion coefficients of twist and writhe respectively, $\theta(x)$ is the Heaviside function, while $\kappa>0$ introduces a surface tension in writhed domains. The reaction terms, called $g$ and $-g$ in the equations for twist and writhe, are equal and opposite so that the linking number, or local catenation number $\sigma(x,t)$ – equal to ${\rm Tw}(x,t)+{\rm Wr}(x,t)$, is globally conserved at all times. We considered $\alpha\geq\beta\geq\gamma\geq 0$, while $\sigma_{0}$ is the local linking/catenation number at which plectonemes appear. Finally, $\sigma_{\eta}$ controls the strength of noise, and $\eta(x,t)$ is a white noise with zero mean and variance equal to

[TABLE]

where $\langle\cdot\rangle$ denotes averaging over different noise realisations. In the following calculations we set $\kappa=0$ for simplicity – this does not affect our analytical treatment.

Let us now assume that Eq. 10 admits a stationary solution, where twist and writhe are constant in steady state. This solution will be such that the reaction term vanishes, $g=0$. Additionally, by summing the two equations we obtain

[TABLE]

Consequently, for a stationary solution, and if there is no twist or writhe flux at the boundaries (i.e., $\frac{\partial{\rm Tw}}{\partial x}=\frac{\partial{\rm Wr}}{\partial x}=0$ there), we obtain that $D_{t}{\rm Tw}+D_{w}{\rm Wr}=C$, where $C$ is a constant which is determined by the initial condition at $t=0$. Intersecting the curve $g=0$ with the straight line $D_{t}{\rm Tw}+D_{w}{\rm Wr}=C$ gives the possible stationary solutions of the system. Notably, for a range of values of $C$ the straight line $D_{t}{\rm Tw}+D_{w}{\rm Wr}=C$ intersects $g=0$ in three points, and the two filled circles in Fig. S8 denotes the two stable points – these are the analog of binodal points for a phase separating system. [Note that we work in the quadrant with ${\rm Tw}>0$ and ${\rm Wr}>0$ without loss of generality, as the rest of the behaviour can be obtained using symmetry.]

The presence of two distinct stable stationary points means that the system will phase separate into two (or more) domains, with a different stable point in each pair of neighbouring domains, and where domain size is set by the global constraint that the total linking number is globally conserved. This argument is the same as that used in Halatek and Frey (2018) to demonstrate the formation of spatially varying stationary patterns in MinD protein systems in bacteria. It is important to note that in order for two stable solutions to exist – equivalently, in order for phase separation to arise – we need $D_{t}>D_{w}$, which is the case for DNA.

It is also instructive to write down an effective equation for the local linking/catenation number, $\sigma={\rm Tw}+{\rm Wr}$, which corresponds to a globally conserved variable. To this end, we can add and subtract the equations for twist and writhe in Eq. 10, to obtain

[TABLE]

where we called $\delta={\rm Tw}-{\rm Wr}$.

Now, we perform a gradient expansion in the second equation to express $\delta$ as a function of $\sigma$. To zero-th order in the gradients, the relation between $\delta$ and $\sigma$ is obtained by setting $g=0$, namely,

[TABLE]

Plugging this into the equation for $\sigma$, we find that it can be written as an effective “model B” dynamics (in the terminology of Chaikin and Lubensky (1995)), which is the relevant equation of motion for a globally conserved order parameter. The equation explicitly reads as follows,

[TABLE]

with $M$ a mobility (which we assume constant), $\mu_{\rm eff}$ an effective chemical potential, and $f_{\rm eff}$ an effective free energy density. The latter is given by:

[TABLE]

In Eq. 16, $A$ is a constant to ensure continuity of $f_{\rm eff}$ at $\sigma={\sigma}_{0}$, whereas given the choices of $\alpha$, $\beta$, $\gamma$ $C_{1}>C_{2}>0$. Therefore the effective free energy is piecewise quadratic, with a singularity at ${\sigma}_{0}$. This function is non-convex so that a common tangent construction predicts phase separation, in line with our previous analysis of the reaction-diffusion model. We note that an effective free energy with the same functional form was proposed by Marko in Marko (2007) to study the thermodynamics of plectoneme formation in supercoiled DNA under tension. Whilst our system of equations leads to phase separation and whilst we can compute the value of the twist and writhe of the coexisting states, a linear stability analysis as in Halatek and Frey (2018) shows that the only point where the system undergoes spinodal decomposition is $\sigma={\rm Tw}+{\rm Wr}={\sigma}_{0}$, which can also be seen in our effective free energy, as $\frac{\partial^{2}f_{\rm eff}}{\partial x^{2}}$ is always positive except at $\sigma={\sigma}_{0}$ – this is due to the choice of a singular reaction term (or equivalently effective free energy). As shown in Fig. 7, the dynamics of the continuum model is similar to that of our molecular dynamics simulations. Additional kymographs are shown in Fig. 9 for larger noise strength. These show some examples of hopping. The mechanism leading to these is similar to that identified in the molecular dynamics simulations, with a plectoneme nucleating approximately at the same time when one of similar size disappears – the combined effect is an apparent hopping of the plectoneme along the chain.