Dynamical Scaling and Phase Coexistence in Topologically-Constrained DNA Melting

TL;DR

This study investigates how the topology of DNA, especially supercoiling, influences its melting behavior, revealing phase coexistence and a topology-dependent scaling law through simulations and theoretical analysis.

Contribution

It introduces a combined simulation and mean field theory approach to explain the broad melting transition and bubble growth in topologically constrained DNA.

Findings

Topologically constrained DNA exhibits phase coexistence during melting.

Melting transition broadens due to competition between melting and supercoiling.

A topology-dependent scaling law governs bubble growth in supercoiled plasmids.

Abstract

There is a long-standing experimental observation that the melting of topologically constrained DNA, such as circular-closed plasmids, is less abrupt than that of linear molecules. This finding points to an intriguing role of topology in the physics of DNA denaturation, which is however poorly understood. Here, we shed light on this issue by combining large-scale Brownian Dynamics simulations with an analytically solvable phenomenological Landau mean field theory. We find that the competition between melting and supercoiling leads to phase coexistence of denatured and intact phases at the single molecule level. This coexistence occurs in a wide temperature range, thereby accounting for the broadening of the transition. Finally, our simulations show an intriguing topology-dependent scaling law governing the growth of denaturation bubbles in supercoiled plasmids, which can be understood within the proposed mean field theory.

Equations

$$\frac{1}{L}{\int }_0^L\sigma (x,t)dx={\sigma }_0\forall t,$$

$$\beta f(\varphi ,\sigma )=$$

$$+a_{\sigma }{\sigma }^2+b_{\sigma }{\sigma }^4+\chi \sigma {\varphi }^2.$$

$$\mu ({\sigma }_{-})=\mu ({\sigma }_{+})$$

$$\Pi ({\sigma }_{-})=\Pi ({\sigma }_{+}).$$

$$\frac{\partial \varphi (x,t)}{\partial t}=-{\Gamma }_{\varphi }\frac{\delta \mathcal{H}(\varphi ,\sigma )}{\delta \varphi }$$

$$\frac{\partial \sigma (x,t)}{\partial t}={\Gamma }_{\sigma }{\nabla }^2\frac{\delta \mathcal{H}(\varphi ,\sigma )}{\delta \sigma }$$

$$\mathcal{H}(\varphi ,\sigma )=\int \left(f(\varphi ,\sigma )+{\gamma }_{\varphi }{\left(\nabla \varphi \right)}^2+{\gamma }_{\sigma }{\left(\nabla \sigma \right)}^2\right)dx,$$

$$\frac{\partial \varphi (x,t)}{\partial t}=$$

$$=-{\Gamma }_{\varphi }\left[2\left(\frac{b^2}{4c}+1-a\right)\varphi +3b{\varphi }^2+4c{\varphi }^3+2\chi \varphi \sigma -{\gamma }_{\varphi }{\nabla }^2\varphi \right]$$

$$\frac{\partial \sigma (x,t)}{\partial t}={\Gamma }_{\sigma }{\nabla }^2\left[2a_{\sigma }\sigma +4b_{\sigma }{\sigma }^3+\chi {\varphi }^2-{\gamma }_{\sigma }{\nabla }^2\sigma \right].$$

$$l(t)\sim \{\begin{array}{l}{t}^1\text{ for topologically free DNA,}\\ t^{1/2}\text{ for topologically constrained DNA.}\end{array}$$

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

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

$$U_{bb}(r)=-\frac{{ϵ}_{bb}{R}_0^2}{2}\log \left[1-{\left(\frac{r}{R_0}\right)}^2\right]\text{ if }r<R_0$$

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

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

$$U_{plan}=\frac{{ϵ}_h}{2}(\alpha -{\alpha }_0{)}^2.$$

$$U_{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 }\mathbf{f}$$

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

$$\beta f(\varphi ,\sigma )=\left(\frac{b^2}{4c}+1-a(T)\right){\varphi }^2+b{\varphi }^3+c{\varphi }^4+a_{\sigma }{\sigma }^2+b_{\sigma }{\sigma }^4+\chi \sigma {\varphi }^2.$$

$$\beta f(\varphi ,\sigma )=\left(3-a(T)\right){\varphi }^2-4{\varphi }^3+2{\varphi }^4+{\sigma }^2+\frac{{\sigma }^4}{2}+\chi \sigma {\varphi }^2.$$

$$a_c(\sigma ,\chi )=\frac{b^2+4c(1+b+c)}{4c}+\chi \sigma$$

$$a_c(\sigma ,\chi )=1+\chi \sigma$$

$${\varphi }_0=\underset{\varphi }{\min }\left\{\frac{\partial f(\varphi ,\sigma )}{\partial \varphi }=0\right\}.$$

$${\varphi }_0(\chi ,\sigma ,a)=\{\begin{array}{l}0\text{ if }a<a_c(\sigma ,\chi )\\ \frac{\left(3+\sqrt{-3+4a-4\chi \sigma }\right)}{4}\text{ if }a\ge a_c(\sigma ,\chi )\end{array}.$$

$${\sigma }_c(a)=\left(a-1\right)/\chi .$$

$${\partial }_{\sigma }^{2}f({\varphi }_0,\sigma )=0.$$

$$a_s(\sigma )=\frac{72{\sigma }^5+27{\sigma }^4+24{\sigma }^3+9{\sigma }^2+2\sigma +3}{36{\sigma }^4+12{\sigma }^2+1}$$

$$\mu (s)={\frac{\partial f({\varphi }_0,\sigma )}{\partial \sigma }|}_s$$

$$\Pi (s)=f({\varphi }_0,s)-{\frac{\partial f({\varphi }_0,\sigma )}{\partial \sigma }|}_{s}s$$

$$\mu ({\sigma }_{-})=\mu ({\sigma }_{+})$$

Full text

One of the most fascinating aspects of DNA is that its biological function is intimately linked to its local topology Bates and Maxwell (2005). For instance, DNA looping Alberts et al. (2014); Rippe et al. (1995) and supercoiling Liu and Wang (1987); Bates and Maxwell (2005); Brackley et al. (2016) are well-known regulators of gene expression, and a variety of proteins, such as Polymerases, Gyrases and Topoisomerases, can affect genomic function by acting on DNA topology Alberts et al. (2014); Bates and Maxwell (2005).

Fundamental biological processes such as DNA transcription and replication are associated with local opening of the double helix, a phenomenon that can be triggered in vitro by varying temperature, pH or salt concentration Wartell and Benight (1985). The melting transition of DNA from one double-stranded (ds) helix to two single-stranded (ss) coils has been intensively studied in the past by means of buoyant densities experiments Vinograd et al. (1968), hyperchromicity spectra Jensen and von Hippel (1976), AFM measurements Jeon et al. (2010), single-molecule experiments Vlijm et al. (2015) and fluorescence microscopy Altan-Bonnet et al. (2003).

In particular, experiments Vinograd et al. (1968); Wartell and Benight (1985) and theories Poland and Scheraga (1966); Kafri et al. (2000) have shown that the “helix-coil” transition in linear or nicked DNA molecules, which do not conserve the total linking number between the two strands, is abrupt and bears the signature of a first-order-like transition. On the other hand, the same transition is much smoother for DNA segments whose linking number is topologically preserved, such as circular, covalently-closed ones Vinograd et al. (1968); Gagua et al. (1981).

Understanding the physical principles underlying DNA melting in topologically constrained (tc) DNA is important since this is the relevant scenario in vivo. For instance, bacterial DNA is circular, while in eukaryotes DNA wraps around histones Alberts et al. (2014), and specialised proteins are able to inhibit the diffusion of torsional stress Naughton et al. (2013).

Intriguingly, and in stark contrast with the behaviour of linear, or topologically free (tf), DNA, the width of the melting transition of tcDNA is relatively insensitive to the precise nucleotide sequence Gagua et al. (1981) thereby suggesting that a universal physical mechanism, rather than biochemical details, may underlie the aforementioned broadening. While biophysical theories of tcDNA melting do exist, they do not reach a consensus as to whether the transition should weaken or disappear altogether Rudnick and Bruinsma (2002); Kabakçioǧlu et al. (2009); it is also unclear if stable coexistence of ss (coiled) and ds (helical) phases may in general be possible regardless of sequence specificity Benham (1979); Sen and Majumdar (1988); Bauer and Benham (1993).

Here we address this topic by employing a combination of complementary methods. First, we perform large-scale coarse-grained Brownian Dynamics (BD) simulations of 1000 base-pairs (bp) long topologically free and constrained double-stranded DNA molecules undergoing melting Fosado et al. (2016). These unprecedentedly large-scale simulations at single bp resolution allow us to measure both kinetic and equilibrium observables of the melting transition. Second, we propose and study a phenomenological Landau mean field theory which couples a critical “denaturation” field ($\phi$) with a non-critical “supercoiling” one ($\sigma$). This approach captures the interplay between local DNA melting and topological constraints, and predicts the emergence of phase coexistence within a wide temperature range, in line with our simulations, and accounting for the experimental broadening of the melting transition in tcDNA. We also derive dynamical equations for the fields $\phi$ and $\sigma$ and discuss the topology-dependent exponents describing the coarsening of denaturation bubbles during DNA melting.

Melting Curves and Phase Diagram –

We first investigate the melting behaviour of DNA by means of BD simulations of the model proposed in Fosado et al. (2016). The dsDNA is made up by two single-stranded chains of “patchy-beads” connected by permanent FENE bonds. Every patch-bead complex represents one nucleotide and complementary strands are paired by bonds connecting patches. We model these bonds as breakable harmonic springs, which mimic hydrogen (H) bonds between nucleotides (see SM for further details). For simplicity, we consider a homopolynucleotide (no sequence dependence) Gagua et al. (1981).

As anticipated, the double-helical structure can be opened up in vitro either by increasing the temperature ($T$), or by increasing pH or salt concentration: both methods effectively reduce the strength of the H bonds, $\epsilon_{\rm HB}$, in between nucleotides. The simulations reported in this Letter emulate the latter route: starting from an equilibrated dsDNA molecule, we perform a sudden quench of $\epsilon_{\rm HB}$, and record the time evolution of the system until a new steady state is reached (see SM Fig. S2).

An observable that directly compares with experiments is the fraction of denatured base-pairs (bp), $\vartheta$. The plot of the equilibrium value $\langle\vartheta\rangle$ as a function of temperature or bond strength can be identified with the melting curve for DNA. Typical profiles obtained from experiments Vinograd et al. (1968) and BD simulations performed in this work, are shown in Fig. 1(A-B): the qualitative agreement is remarkable.

The sharpness of the melting transition can be quantified in terms of the maximum value attained by the differential melting curve as $\Delta t=\left|d\langle\vartheta\rangle/dt\right|^{-1}$, where $t$ can either be temperature, $T$, or effective hydrogen bond strength, $\epsilon_{\rm HB}$, depending on the denaturation protocol. Quantitatively, Figures 1(A)-(B) show that experiments and simulations agree in predicting melting curves for tcDNA about three times broader than for tfDNA, i.e. $\left.\Delta t\right|_{\rm tc}/\left.\Delta t\right|_{\rm tf}\simeq 3$.

From these observations, it is clear that the melting behaviour of DNA is affected by global topology. On the other hand, melting occurs through local opening of the double-helical structure. The challenge faced by a theory aiming to understand the “helix-coil” transition in tcDNA is therefore to capture local effects due to the global topological invariance. To this end, it is useful to define an effective local supercoiling field $\sigma(x,t)\equiv\left(Lk(x,t)-Lk_{0}\right)/Lk_{0}$, where $Lk_{0}$ is the linking number between the two strands in the relaxed B-form state, i.e. 1 every 10.4 bp, and $Lk(x,t)$ is the effective linking number at position $x$ and time $t$. For a circular closed molecule of length $L$

[TABLE]

where $\sigma_{0}$ is the initial supercoiling deficit, which can be introduced and locked into the chain by, for instance, the action of topological enzymes Alberts et al. (2014). On the contrary, circular nicked or linear (tf) dsDNA molecules need not satisfy Eq. (1), since any deviation from the relaxed supercoiling state can be expelled through the chain ends or the nick. In light of this, it is clear that subjecting a tcDNA to denaturation-promoting factors causes a competition between entropy and torsional stress: the former associated with the denatured coiled regions Poland and Scheraga (1966); Kafri et al. (2000), the latter arising in the intact helical segments Kabakçioǧlu et al. (2009).

Motivated by these observations, we propose the following phenomenological mean field theory for the melting of tf and tcDNA. We consider a denaturation field, $\phi(x,t)$, describing the state of base-pair $x$ at time $t$ (e.g., taking the value [math] if intact or $>0$ if denatured), coupled to a conserved field, $\sigma(x,t)$, tracking the local supercoiling. A Landau free energy can be constructed by noticing that: (i) the denaturation field $\phi$ should undergo a first-order phase transition when decoupled from $\sigma$ Poland and Scheraga (1966), (ii) the elastic response to the torsional stress should be associated with even powers of $\sigma$ Benham (1979) 111This is a good approximation for small $\sigma$., and (iii) the coupling between $\sigma$ and $\phi$ should be such that there should be an intact dsDNA phase at sufficiently low $T$, i.e. $\phi=0$ for any $\sigma_{0}$ at $T<T_{c}$.

Based on these considerations we can write an effective free energy density as:

[TABLE]

In Eq. (Melting Curves and Phase Diagram –), the first term is written so that the parameter $a(T)\sim T/T_{c}$, $b<0$ and we keep the quartic term in $\sigma$ to ensure there is a local minimum around $\sigma=-1$ (or $Lk(x,t)=0$) when $\phi=\phi_{0}$, corresponding to the fully denatured state (see SM). Finally, the coupling term $\chi\sigma\phi^{2}$ models the interplay between supercoiling and local melting; it can be turned off by simply setting $\chi=0$ to approximate tfDNA (in this framework the torsional stress can be expelled infinitely fast). We also highlight that by setting $\chi=0$, Eq. (Melting Curves and Phase Diagram –) predicts a first-order melting transition.

The free energy density in Eq. (Melting Curves and Phase Diagram –) displays a minimum at $\phi=0=\phi_{\rm ds}$ (helical state), and can develop a competing one at $\phi=\phi_{\rm ss}>1$ (coiled state) which in general depends on $\sigma$, $\chi$ and $a$ (see SM). The free energy density of the two becomes equal at the critical temperature $a=a_{c}$, which reduces to $a_{c}(\sigma)=1+\chi\sigma$, by fixing $c=-b/2$ (see SM). This relation states that more negatively supercoiled molecules denature at lower temperatures, as in experiments Víglaský et al. (2000), whereas for tfDNA ($\chi=0$) the critical temperature is independent on supercoiling.

To obtain the phase diagram of the system in the space $(a,\sigma_{0})$ we focus on dsDNA molecules with fixed, and initially uniform, value of supercoiling $\sigma_{0}$, at fixed temperature $a$. For such conditions, the system attains its free-energy minimising state for a value of $\phi=\phi_{0}(\chi,\sigma_{0},a)$ Matsuyama et al. (2002). The uniform solution $(\sigma_{0},\phi_{0})$ is linearly unstable if it lies within the spinodal region (in Fig. 2 shaded in grey), i.e. where $\partial^{2}f(\phi_{0},\sigma)/\partial\sigma^{2}\leq 0$ Chaikin and Lubensky (2000). In Figure 2 we fix for concreteness $b=-4$, $a_{\sigma}=1$, $b_{\sigma}=1/2$, $\chi=2$ (different parameter choices lead to similar diagrams provided $b$ remains negative 222Given these values of $b$, $a_{\sigma}$, $b_{\sigma}$, we further require $1.72\leq\chi\leq 2.67$ to obtain a spinodal region.).

A system with unstable uniform solution separates into two phases with low ($\sigma_{-}$) and high ($\sigma_{+}$) supercoiling levels, as this lowers the overall free energy. The values of $\sigma_{-}(a)$ and $\sigma_{+}(a)$ are the coexistence curves, or binodals, which are found by imposing that both chemical potential $\mu(s)\equiv\left.\partial f(\phi_{0},\sigma)/\partial\sigma\right|_{s}$ and pressure $\Pi(s)=f(\phi_{0},s)-\mu(s)s$ must be equal in the two phases Matsuyama et al. (2002). This translates into solving a system of two equations with two unknowns,

[TABLE]

By noticing that $\sigma_{0}$ needs to satisfy Eq. (1) for tcDNA, it is straightforward to find the fractions of the system in the high and low supercoiling phases as $f_{+}=(\sigma_{0}-\sigma_{-})/(\sigma_{+}-\sigma_{-})$ and $f_{-}=(\sigma_{+}-\sigma_{0})/(\sigma_{+}-\sigma_{-})$, respectively.

The phase diagram in the $(a,\sigma_{0})$ space is reported in Figure 2, where we show that the coexistence lines $\sigma_{-}(a)$ and $\sigma_{+}(a)$ wrap around the critical first-order transition line $\sigma_{c}(a)=(a-1)/\chi$ which therefore becomes “hidden” Matsuyama et al. (2002). In light of this we argue that the smoother transition observed for tcDNA Vinograd et al. (1968); Gagua et al. (1981) can be understood as a consequence of the emergence of a coexistence region in the phase space which blurs the underlying first-order transition. This argument also explains the “early melting” of closed circular DNA Vinograd et al. (1968), which can be understood as the entry into the coexistence region from low temperatures.

Intriguingly, our phase diagram includes a region (cross-shadowed in Fig. 2) where the system displays stable coexistence of two open phases (i.e. $\phi=\phi_{\rm ss}$ in both sub-systems) with supercoiling levels $\tilde{\sigma}_{-}$ and $\tilde{\sigma}_{+}$.

Dynamical Scaling –

The dynamics of the non-conserved order parameter, $\phi$, and the conserved one, $\sigma$, can be found following the Glauber and Cahn-Hilliard prescriptions, respectively. Consequently, the system can be described by the following “model C” equations Chaikin and Lubensky (2000)

[TABLE]

where $\Gamma_{\sigma,\phi}$ are relaxation constants,

[TABLE]

and $\gamma_{\phi,\sigma}$ determine the effective surface tension of bubbles and supercoiling domains, respectively. From Eqs. (Melting Curves and Phase Diagram –) and (5) one can write

[TABLE]

We numerically solve this set of partial differential equations (PDE) on a 1D lattice of size $L$ for fixed $a$ and $\sigma_{0}$ (see SM for details) and compare the evolution of denaturation bubbles with the one observed in BD simulations. Note this set of equations disregards thermal noise, hence it is in practice a mean field theory.

In Figure 3 we show “kymographs” from BD simulations, capturing the state of each base-pair (either intact or denatured) against time for tf and tcDNA. As one can notice, after the energy quench at $t=0$, the linear (tfDNA) molecule starts to denature from the ends and eventually fully melts. On the other hand, in the closed circular (tcDNA) molecule, bubbles pop up randomly over the whole contour length, and the steady state entails a stable fraction, $0<\vartheta<1$, of denatured bp (see also SM, Fig. S2).

We observe a similar behaviour when the fields $\phi$ and $\sigma$ are evolved via Eqs. (7), starting from a single small bubble at temperature $a$ within the coexistence region (see Fig. 3D and SM). While the bubble grows, the supercoiling field is forced outside the denatured regions and accumulates in the ds segments. The increasing positive supercoiling in the helical domains slows down and finally arrests denaturation, resulting in phase coexistence in steady state, between a denatured phase with $\sigma=\sigma_{-}$ and $\phi=\phi_{\rm ss}>1$, and an intact phase with $\sigma=\sigma_{+}$ and $\phi=\phi_{\rm ds}=0$.

The growth, or coarsening, of a denaturation bubble, $l$, can be quantified within our mean field theory and BD simulations: in the former case by numerical integration of Eq. (7), in the latter by measuring the size of the largest bubble over time and averaging over independent realisations. As shown in Fig. 4(A-C) we find that in both mean field and BD simulations,

[TABLE]

In other words, we find that the exponent $\alpha$ governing the local growth of a denaturation bubble depends on the global topology of the molecule.

We propose the following argument to explain the values of $\alpha$. For tfDNA (e.g., nicked or linear), we can assume that the supercoiling field relaxes quickly, and gets expelled outside, without affecting the dynamics of the denaturation field. In this case, the free energy can be approximated as $f\simeq(\epsilon_{\rm HB}-T\Delta S)l$, so that there is a constant increase in entropy per each denatured bp when $T>T_{c}=\epsilon_{\rm HB}/\Delta S$. This implies that Chaikin and Lubensky (2000) $\psi dl/dt\simeq df/dl\sim const$, with $\psi$ an effective constant friction; as a result we obtain $l(t)\sim t$.

On the other hand, the value of $\alpha=1/2$ observed for tcDNA (e.g., circular non-nicked plasmids) can be understood by quantifying the slowing down of denaturation due to the accumulation of a “wave” of supercoiling, raked up on either side of the growing bubble. We argue that the flux of $\phi$ through a base pair at the bubble/helix interface is $J_{\phi}\sim\phi_{\rm ss}dl/dt$. At the same time, the flux of $\sigma$ can be obtained by noticing that the “wave” can be approximated by a triangle with constant height $h=\sigma_{+}-\sigma_{0}$ and base $b\sim l(t)$ (see SM Fig. S10 and Fig. 3(C)). This is because the total supercoiling enclosed by the wave must be proportional to the one expelled from within the denatured bubble, which is $\sim\left|\sigma_{-}-\sigma_{0}\right|l(t)$. One can therefore write that $J_{\sigma}=-\Gamma_{\sigma}\partial_{x}\sigma\simeq\Gamma_{\sigma}h/l$, for the supercoiling flux in, say, the forward direction. At equilibrium, the two fluxes must balance, i.e. $J_{\phi}\sim J_{\sigma}$, and therefore $\phi_{\rm ss}dl/dt\sim\Gamma_{\sigma}h/l$, or $l(t)\sim t^{1/2}$. Finally, we highlight that this argument depends on the slowing down over time of 1D supercoiling fluxes, hence is qualitatively distinct from the reason why $\alpha=1/2$ in dimensions $d\geq 2$ in (non-conserved) model A Chaikin and Lubensky (2000).

Linking within Denaturation Bubbles –

As a final result, we perform BD simulations to characterise the topology of a denaturation bubble through the linking number that can be stored inside it (Fig. 4(D)). An idealised bubble is identified by ${\rm Lk}=0$ (and $\sigma=-1$) Kabakçioǧlu et al. (2009), whereas our BD simulations show that a denatured region of fixed size $l$ (imposed by selectively breaking only $l$ consecutive bonds along an intact dsDNA molecule) has a non-zero linking number $Lk_{d}$ (see Fig. 4(D)) 333In practice, $Lk_{d}$ was obtained by removing the linking number of the chain outside the denatured region $Lk_{\rm out}$ from the initially set $Lk$, see SM for details..

We found that for small $l$, $Lk_{d}$ displays a remarkable signature of global topology (through the value of $\sigma_{0}$); instead, the scaling behaviour at large $l$ appears to follow $Lk_{d}\sim l^{1.25}$ irrespectively of $\sigma_{0}$, until it reaches $Lk_{0}$. The finding that a denaturation bubble displays a non-trivial and $l$-dependent linking number suggests that idealised ($Lk_{d}=0$) bubbles may not always be reflecting realistic behaviour. Further, it may be of relevance for processes such as DNA replication, as it suggests that supercoiling or torsional stress may be able to diffuse past branching points such as replication forks Postow et al. (2001).

Conclusions –

In summary, we have studied the melting behaviour of topologically constrained DNA through a combination of large-scale BD simulations and mean field theory. A key result is that the phase diagram for tcDNA melting generally involves a phase coexistence region between a denatured and an intact phases, pre-empting a first-order denaturation transition as in tfDNA. This finding provides a theoretical framework to explain the long-standing experimental observation that the denaturation transition in circular, and not nicked, supercoiled plasmids is seemingly less cooperative (smoother) than for linear, or nicked, DNA Vinograd et al. (1968); Gagua et al. (1981).

We have further studied, for the first time, the coarsening dynamics of denaturation bubbles in tcDNA, and found a remarkable agreement between BD simulations and mean field theory, both reproducing similar topology-dependent scaling exponents that can be understood within our theoretical model. It would be of interest to investigate such dynamics experimentally in the future.

DMa and DMi acknowledge ERC for funding (Consolidator Grant THREEDCELLPHYSICS, Ref. 648050). YAGF acknowledges support form CONACyT PhD grant 384582.

I Brownian Dynamics Simulations of double stranded DNA

Brownian Dynamics simulations are performed using the model developed in Ref. Fosado et al. (2016) to which we refer for details on parameters and validation. Below we briefly describe its main features.

The double-stranded (ds) DNA molecule is modelled as two polymers made of “patchy-beads” connected by springs. Specifically, each nucleotide is represented by a rigid body made up of two spherical monomers: a bead which models the sugar-phosphate backbone, and a “patch” that represents the nitrogenous base. Two nucleotides belonging to complementary strands are shown in Figure S1(a). The beads have an excluded volume of $\sigma=1$ nm in diameter and they are shown in red for one strand and in blue for the complementary. The corresponding patches are shown in cyan and pink: these have no associated excluded volume and are placed at a distance of $\sigma/2$ from the centre of the bead. Excluded volume between beads is modelled via a truncated and shifted Lennard-Jones (LJ) potential (also known as Weeks-Chandler-Anderson potential),

[TABLE]

if $r<2^{1/6}\sigma$, and $U_{LJ}(r)=0$ otherwise.

The hydrogen bonding interaction, responsible for holding two nucleotides belonging to complementary strands together, is modelled by breakable harmonic springs acting between two patches (Fig. S1(a)). The potential associated to this interaction is

[TABLE]

if $r\leq r_{\rm c,hb}$ and $U_{\rm hb}=0$ otherwise. The equilibrium bond distance $r_{\rm 0,hb}$ is set to zero and the critical distance $r_{\rm c,hb}=0.3$ nm marks the point at which the bond breaks. The strength of the potential is $\epsilon_{HB}$, the main parameter used in the main text to drive the melting transition.

A single strand is formed by a chain of patchy-beads connected via FENE bonds of length $d_{bp}=0.46$ nm. This potential reads

[TABLE]

with $R_{0}=0.68\sigma$ and $U_{bb}(r)=\infty$ otherwise.

We also model the stacking of base-pair by means of a Morse potential which keeps the distance between consecutive patches along each strand around $r_{\rm 0,stack}=0.34$ nm, i.e.

[TABLE]

While the choice of the two previous potentials ensures that the pitch of the chain is the one expected for dsDNA (around 10 bp per helical turn), a dihedral interaction (see Fig. S1(b)) between the particles of two consecutive nucleotides in the same strand enforces right-handed helicity for the DNA molecule – i.e., we set

[TABLE]

where $d=180^{\circ}-36^{\circ}$.

Planarity between basis is imposed through a harmonic potential preventing the angle formed between two consecutive patches and one bead (all in the same strand) to change from its equilibrium value ($\alpha_{0}=90^{\circ}$).

[TABLE]

Finally, in order to regulate the stiffness of the chain we use a Kratky-Porod potential controlling the angle between three consecutive patches along one strand, i.e.

[TABLE]

A full description of all the potentials used in the model can be found in Ref. Fosado et al. (2016).

In order to avoid large overlaps between beads and to preserve the correct geometry of DNA we consider two types of beads in each strand: sterically interacting beads (shown as small solid red spheres for one strand and blue for the other in Fig. S1(C)) are intercalated by two non-sterically interacting beads (represented as small grey spheres). These do not interact with the beads along the same strand but they repel all beads on the complementary strand within an excluded radius $\sigma/2$. This choice ensures that only non-overlapping beads along each single strand sterically interact with one another and, at the same time, allows us to preserve the topology by forbidding the two strands crossing through one another.

The nucleotides made up by the bead-patch complex evolve as a rigid bodies undergoing Langevin dynamics, i.e. the position of each bead $\bm{r}$ is updated according to equation

[TABLE]

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

[TABLE]

along each Cartesian component (denoted by Greek letters). Integration of Eq. (8) is performed using a velocity-Verlet algorithm through the LAMMPS Plimpton (1995) code, run in BD mode.

II Dynamics of Melting

In the main text we show the melting profiles $\langle\vartheta\rangle$ as a function of the effective bond strength $\epsilon_{\rm HB}$ in Figure 1. The value of stably denatured bp $\langle\vartheta\rangle$ is obtained by averaging the last $2$ $10^{7}$ $\tau_{LJ}$ timesteps of the trajectories obtained for $\vartheta(t)$ (examples shown in Fig. 2). As one can notice from the figure, all the samples start from non-denatured states $\vartheta(0)=0$; after a sudden quench in $\epsilon_{\rm HB}$ the systems evolve until a new steady state is reached. For topologically constrained (closed circular) molecules, $\vartheta(t)$ reaches a steady state value in between 0 and 1 (see red, green and blue curves). On the other hand, linear molecules initially denature very little but eventually (after a time 10 times longer than the equilibration time for rings) fully denature.

A

III Landau Free-Energy

Here we study the free energy whose density is written in Eq. (Melting Curves and Phase Diagram –),

[TABLE]

In this equation, $\phi=\phi(x,t)$ is the field describing the state (coiled or helical) of base-pair $x$ along the contour of the DNA molecule at time $t$ and $\sigma=\sigma(x,t)$ the field quantifying the amount of local supercoiling. The parameter $a=a(T)\sim T/T_{c}$ is the only temperature-dependent parameter and $\chi$ the coupling between $\phi$ and $\sigma$.

This free energy density describes a first-order transition between a closed (double-stranded, ds, $\phi=\phi_{\rm ds}=0$) and an open (single stranded, ss, $\phi=\phi_{\rm ss}>1$) DNA molecule. We choose the parameter $c=-b/2$ so that the critical temperature depends only on $\chi\sigma$ (see below) and so that the free energy density displays a minimum located at around $\sigma=-1$ when $\phi$ attains its free-energy minimising value of $\phi_{0}$ (see below). The corresponding choices are $c=-b/2$, $a_{\sigma}=1$, $b_{\sigma}=1/2$ and $b=-4$. The free energy density then reduces to

[TABLE]

As a function of $\phi$, Eq. 11 can display two local minima whose energy is equal for topologically unconstrained ($\chi=0$) molecules at the critical (melting) temperature $T_{c}$ (defined via $a(T_{c})=a_{c}=1$). In this situation (see Fig. S3) the minima are located at $\phi_{\rm ds}=0$ (closed, double stranded) and at $\phi_{\rm ss}=1$ (open, single stranded) and mark the equal probability coexistence point.

In general, the critical temperature is supercoiling-dependent through the coupling parameter $\chi$ (see Fig. S4A-C). The function $a_{c}(\sigma,\chi)$ is simply found by solving $f(\phi_{\rm ds},\sigma)=f(\phi_{\rm ss},\sigma)$, i.e.

[TABLE]

or

[TABLE]

by substituting the values of parameters. We report $a_{c}$ in Fig. S4(D) for $\chi=0$ and $\chi=2$. It can be readily appreciated that for topologically constrained ($\chi>0$) molecules, the critical temperature $a_{c}$ depends linearly on the supercoiling, in agreement with experimental observations Víglaský et al. (2000).

The equilibrium value of the denaturation field $\phi_{0}$ is a discontinuous function of the temperature and is found by minimisation of the free energy, as

[TABLE]

This condition leads to (see also curves in Fig. S4(C))

[TABLE]

Because $\phi$ is a non-conserved order parameter, it will always attain the value $\phi_{0}(\chi,\sigma,a)$ for any given set of $\chi$, $\sigma$ and $a$. The free energy density of the system can therefore be written as $f(\phi_{0}(\chi,\sigma,a),\sigma)=f(\phi_{0},\sigma)$. The behaviour of this function is plotted in Figure S5.

One can readily check that for $\chi=0$ the discontinuity is always located at $a_{c}=1$, irrespective of $\sigma$, and that the profile of the free energy is symmetric with respect to $\sigma$. On the contrary, for topologically constrained ($\chi>0$) molecules, the critical temperature is supercoiling-dependent, $a_{c}=a_{c}(\sigma)$, and the profile as a function of $\sigma$ is more complex. The term $\sigma^{4}$ in the free energy and the parameters chosen ensure that $f(\phi_{0},\sigma)$ displays a minimum around $\sigma=-1$ for $\phi_{0}=\phi_{\rm ss}$ which corresponds to a fully denatured state with no supercoiling ($Lk=0$) within the coiled region.

Figure S5(D) can be understood as follows: for fixed temperature $a$ the critical supercoiling $\sigma_{c}$ can be found through eq. (13) as

[TABLE]

At $\sigma_{c}(a)$ the system undergoes a first-order transition (corresponding to the discontinuity of the curves in Fig. S5(D)), and the system switches equilibrium state from $\phi_{\rm ds}$ to $\phi_{\rm ss}$. This corresponds to a cusp in the behaviour of $f(\phi_{0},\sigma)$. In order to find the coexisting phases we now need to perform the common tangent construction on the function $f(\phi_{0},\sigma)$ represented in Figure S5(D)).

A

A

IV Spinodal Region and Binodal Lines

The supercoiling field $\sigma$ is a globally conserved field, and the system may therefore undergo phase decomposition in regions of high and low $\sigma$ in order to lower the overall free energy while preserving the initially set $\sigma=\sigma_{0}$. The region in the parameter space $(a,\sigma_{0})$ for which this occurs is know as coexistence region.

The region where the uniform phase $\sigma(x,t)=\sigma(t)$ is linearly unstable is identified as the “spinodal” region and can be readily found by identifying the location of the inflexion points of $f(\phi_{0},\sigma)$ (see Fig. S6). In other words one needs to solve

[TABLE]

This condition gives an analytical expression for the spinodal lines (for $\chi=2$) as (see Fig.S 7)

[TABLE]

valid if $-0.57735<\sigma_{0}<0.57735$.

The range of values of $\chi$ for which the requirement in Eq. (17) returns a solution $a_{s}(\sigma)$ is $1.72\leq\chi\leq 2.67$. For this reason we chose a value in this range, i.e. $\chi=2$.

When the system undergoes phase decomposition in the conserved field $\sigma$, two coexisting solutions (or phases) appear, with distinct values of supercoiling: we label $\sigma_{+}$ the higher supercoiling solution and $\sigma_{-}$ the lower supercoiling transition (in general $\sigma_{+}$ and $\sigma_{-}$ need not have opposite sign). he values of $\sigma_{+}$ and $\sigma_{-}$ are found by the common tangent construction, i.e. by equating the chemical potentials

[TABLE]

and the pressures $\Pi(s)$

[TABLE]

in the two phases, thereby leading to a system of two equations in two unknowns:

[TABLE]

Equation (21) states that the chemical potentials of the two phases must balance in order for the system to be in equilibrium. Equation (22) restricts the translational degree of freedom and it identifies the only pair of points of $f(\phi_{0},\sigma)$ that have equal tangent $\mu(\sigma_{+})=\mu(\sigma_{-})$ and that can be joined by a straight line.

Graphically, this procedure can be summarised as in Figure S8. The blue curve is the free energy $f(\phi_{0},\sigma)$ resulting from a choice of $a=1.1$. The green line corresponds to $\mu(\sigma_{-})\sigma=\mu(\sigma_{+})\sigma$ and the orange curve to the pressure $\Pi(\sigma)$. The common tangent construction identifies the points $\sigma_{-}$ and $\sigma_{+}$ with same tangent $\mu$ and with equal pressure, i.e. where the value of the tilted free energy density (orange curve) is the same. The red dot placed along the (blue) curve $f(\phi_{0},\sigma)$ represents the value of the free energy density $f(\phi_{0},\sigma_{0})$ for $\sigma_{0}=-0.1$, which is higher than the free energy per unit volume of the phase-separated system:

[TABLE]

represented by the black dot. In eq. (23), $f_{+}$ and $f_{-}$ are the fractions of the system in the high and low supercoil phases, respectively.

For a general value of the initial supercoiling field $\sigma_{0}$, one can find $f_{+}$ and $f_{-}$ through the so-called lever-rule:

[TABLE]

which gives

[TABLE]

and

[TABLE]

Equations (21) and (22) can be solved numerically and give the binodal lines $\sigma_{+}(a)$ and $\sigma_{-}(a)$. These are plotted in Fig. S9 along with the critical line $\sigma_{c}(a)$ and the spinodal line from eq. (18).

As one can notice, the binodal line $\sigma_{-}(a)$ crosses through the spinodal region leaving a reentrance where the uniform solution is linearly unstable, yet the binodals $\sigma_{+}(a)$ and $\sigma_{-}(a)$ are not the stable phases. In this region two stable open (ssDNA) phases are observed, whose supercoiling take values on additional binodals $\tilde{\sigma}_{+}$ and $\tilde{\sigma}_{+}$ shown in dashed lines.

V Equations of Motion: Model C

The equations describing the dynamics of the non-conserved field $\phi$ and the conserved one $\sigma$ can be derived from the Landau free energy as Chaikin and Lubensky (2000)

[TABLE]

where $\Gamma_{\phi}$ and $\Gamma_{\sigma}$ are relaxation constant, and $\mathcal{H}$ denotes the total free energy,

[TABLE]

The equations of motion are therefore:

[TABLE]

Equations (30) are solved numerically on a 1D lattice of size $L=1000$ via Euler method with integration time $dt=10^{-3}$ $t$ and lattice spacing $dl=1$. The parameters are set as before: $b=-4$, $c=-b/2$, $a_{\sigma}=1$ and $b_{\sigma}=1/2$. The surface tensions $\gamma_{\sigma}$ and $\gamma_{\phi}$ are both set to $10$ in units of squared lattice spacing. The mobility $\Gamma_{\sigma}$ is set so that the supercoiling field $\sigma$ relaxes more quickly than $\phi$ (as we expect to be the case for DNA during melting), i.e. $\Gamma_{\sigma}=10$ (units $dl^{2}t^{-1}$) while $\Gamma_{\phi}=1$ (units $t^{-1}$).

The system is initialised with a small denatured bubble of size $l_{b}=10\ll L$ so that

[TABLE]

while the supercoiling field is set to

[TABLE]

Finally, since the molecule is circular-closed, we set periodic boundary conditions along $L$ for the values of the fields and their derivatives, i.e. $\phi(-L/2,t)=\phi(L/2,t)$, $\partial_{x}\phi(-L/2,t)=\partial_{x}\phi(L/2,t)$ and $\sigma(-L/2,t)=\sigma(L/2,t)$, $\partial_{x}\sigma(-L/2,t)=\partial_{x}\sigma(L/2,t)$ for any time $t$.

Importantly, we can model topologically unconstrained molecules by setting $\chi=0$, for which the denaturation field evolves uncoupled from $\sigma$. This scenario can be viewed as an approximation for nicked molecules where the supercoiling relaxes infinitely fast (through strand rotations at the nick).

Starting the simulations from the above-mentioned conditions, we can evolve the system and record the dynamics of the fields in time. The typical profiles are reported in Fig. 10(A,B), where the denaturation field $\phi$ and supercoiling field $\sigma$ are shown for $a=1.5$, $\chi=2$ and $\sigma_{0}=-0.06$. For these values of $\sigma_{0}$ and $a$, the uniform solution is (linearly) unstable and for this reason the initially small denatured bubble grows in time and eventually arrests at the point where two stable phases coexist. At the same time, the supercoiling field is “emptied” from inside the bubble, where it takes the stable value of $\sigma_{-}$, and “poured” in the $L-l(t)$ portion of the chain, where it eventually attains the value $\sigma_{+}$.

As discussed in the main text, the size of the denatured bubble $l(t)$ has been found to scale in time as

[TABLE]

The exponent for unconstrained DNA, such as nicked or linear, can be simply explained by assuming that for temperatures above a critical $T_{c}=\epsilon/\Delta S$ the system can lower its (free) energy for each un-paired (denatured) base-pair as

[TABLE]

This directly implies that Chaikin and Lubensky (2000)

[TABLE]

and therefore $l(t)\sim t$.

On the other hand, the value of $\alpha$ observed for topologically constrained, such as circular closed, molecules can be understood by quantifying the slowing down of the denaturation field due to the accumulation of a “wave” of supercoiling on either side of the growing region. In particular, one can argue that the flux of $\phi$ through a base pair at the interface of a growing bubble is

[TABLE]

while the flux of $\sigma$ can be obtained by noticing that the “wave” roughly assumes a triangular shape with height $h=\sigma_{+}-\sigma_{0}$ and base $b\sim l(t)$, so that the area enclosed is proportional to the one expelled from inside the denatured bubble $\sim l(t)$ (see Fig. S10)). In light of this one can approximate

[TABLE]

for, e.g., the flux in the forward direction. At equilibrium, the two fluxes must balance, i.e. $J_{\phi}\sim J_{\sigma}$, and therefore

[TABLE]

or $l(t)\sim t^{1/2}$.

Because the approximation $\left.\partial_{x}\sigma(x,t)\right|_{\rm front}\simeq(\sigma_{+}-\sigma_{0})/l(t)$ might look crude, we numerically computed the minimum and the average slope attained by the $\sigma$ field on the front of the supercoiling wave. These two quantities are plotted in Fig. S11 against the size of the bubble $l$. From the figure one can readily notice that the approximation is in fact not far from the real behaviour.

VI Computing Linking number in denatured bubbles

For a topologically constrained dsDNA the linking number $Lk$ has to be preserved at any time. In particular, even when a circular closed molecule of DNA is fully denatured, all the initial linking number has still to be preserved. One may imagine this situation as an extreme case where the only denaturation bubble takes up the whole DNA molecule. In light of this, one may ask whether a generic bubble of size $l$ can store some non-zero linking number, $Lk_{d}$. This case is generally discarded in simple models Kabakçioǧlu et al. (2009); Jeon et al. (2010); Poland and Scheraga (1966), which consider no supercoiling within denatured regions. On the contrary, arguments from elasticity theory Benham (1979) suggest and support the idea that some linking may still be present within denatured bubbles, since a complete expulsion of this would then contribute to increase the torsional (twist) or bending (writhe) energies of the intact (non-denatured) segments through the well-known $Lk=Wr+Tw$ relation Fuller (1978).

In order to test this argument we performed simulation in which we kept the strength of the hydrogen bond $\epsilon_{HB}=6.0k_{B}T$ constant and high enough so denaturation could not take place. Then, for rings with different level of supercoiling ($\sigma=0.06,0,-0.06$, or $Lk=106,100,94$, respectively), we created a bubble of size $l$ by manually deleting the corresponding consecutive hydrogen-bonds along the chain.

After equilibration we computed the linking number by numerically calculating

[TABLE]

between curves $C_{R}$ and $C_{B}$. In Eq. (40), $\bm{r}_{R}$ and $\bm{r}_{B}$ are the vectors defining the position of the beads along the curves $C_{R}$ and $C_{B}$, respectively.

In order to find the linking number inside the denatured region, $Lk_{d}$, we first replaced the bubble by straight lines joining the first non-denatured base-pairs (see Fig. S12), so that Eq. (40) returned the linking number outside the denatured region, $Lk_{\rm out}$. Finally, we simply use the invariance of the total linking number to find $Lk_{d}=Lk-Lk_{\rm out}$ (see Fig. S13).