A polymer model of bacterial supercoiled DNA including structural transitions of the double helix

TL;DR

This paper develops a comprehensive polymer model of bacterial supercoiled DNA that incorporates structural transitions, enabling better prediction and analysis of DNA conformations and thermodynamics under various forces.

Contribution

It introduces an extended self-avoiding rod-like chain model that accounts for DNA structural variations, validated by simulations and analytical approximations.

Findings

Model accurately captures multiscale properties of long DNA chains.

High-force estimations of DNA parameters can be misleading without proper analysis.

Reevaluation of existing data questions previous parameter assumptions for DNA forms.

Abstract

DNA supercoiling, the under or overwinding of DNA, is a key physical mechanism both participating to compaction of bacterial genomes and making genomic sequences adopt various structural forms. DNA supercoiling may lead to the formation of braided superstructures (plectonemes), or it may locally destabilize canonical B-DNA to generate denaturation bubbles, left-handed Z-DNA and other functional alternative forms. Prediction of the relative fraction of these structures has been limited because of a lack of predictive polymer models that can capture the multiscale properties of long DNA molecules. In this work, we address this issue by extending the self-avoiding rod-like chain model of DNA so that every site of the chain is allocated with an additional structural degree of freedom reflecting variations of DNA forms. Efficient simulations of the model reveal its relevancy to capture multiscale properties of long chains (here up to 21 kb) as reported in magnetic tweezers experiments. Well-controlled approximations further lead to accurate analytical estimations of thermodynamic properties in the high force regime, providing, in combination with experiments, a simple, yet powerful framework to infer physical parameters describing alternative forms. In this regard, using simulated data, we find that extension curves at forces above 2 pN may lead, alone, to erroneous parameter estimations as a consequence of an underdetermination problem. We thus revisit published data in light of these findings and discuss the relevancy of previously proposed sets of parameters for both denatured and left-handed DNA forms. Altogether, our work paves the way for a scalable quantitative model of bacterial DNA.

Figures

Equations

$$E_{\text{sRLC}}(\mathcal{C})=E_B(\mathcal{C})+E_T(\mathcal{C})-fz,$$

$$E_{2sRLC}(\mathcal{C})=E_{\gamma }(\mathcal{C})+E_B(\mathcal{C})+E_T(\mathcal{C})+E_J(\mathcal{C})-fz,$$

$$E_{\gamma }^{eff}(\mathcal{C})$$

$$E_T^{eff}(\mathcal{C})$$

$$F_{f,\sigma }(\lambda )$$

Full text

††journal: Physica A

1 Introduction

In bacteria, and probably in numerous eukaryotes, genome-wide coordination of transcription primarily relies on the physics of DNA in interaction with RNA polymerases [1, 2]. Understanding bacterial transcription coordination therefore requires understanding the physics of DNA, more specifically the physics of under and over-wound DNA [3, 4, 5, 6, 7, 8], also known as DNA supercoiling. DNA is indeed continually processed in cells by topoisomerases [9], whose activity allows relaxing the transient constraints generated by replication and by transcription itself, both processes tending to overwind DNA downstream and to underwind it upstream [10, 11]. As a result, bacterial DNA is generally underwound, which is commonly referred to as the negative supercoiling of bacterial genomes. In model organisms, and probably in most bacteria [12], this negative supercoiling is further partitioned into approximately 10 kb long topologically independent domains [13, 14].

Although supercoiling is produced at specific sites (e.g. at the target sites of gyrases), it may impact DNA globally as a consequence of the conservation of the linking number. This linking number, Lk, is equal to the sum of the twist (Tw), the cumulative helicity of the molecule, plus the writhe (Wr), the global intricacy of the molecule [15]. Thus, for topologically constrained DNA as in the case of circular molecules (e.g. plasmids), of topologically constrained linear domains [10], or when DNA molecules are manipulated by magnetic tweezers [16, 17, 18], any local variation of the twist results in a global variation of the writhe, and reciprocally. As a consequence, negative DNA supercoiling leads to both superstructuring participating to the global compaction of bacterial chromosomes (through e.g. the formation of so-called plectonemes) and to local modifications of the double helix structure. The latter can induce the formation of functionally important DNA forms different from B-DNA [19, 20], including denaturation bubbles [21], cruciforms [22], left-handed Z-DNA [23] and left-handed L-DNA [24, 25, 26, 27].

One of the most important challenges raised by the biophysical and functional characterizations of bacterial genomes thus consists in predicting the multiscale distribution of supercoiling constraints. Such prediction nevertheless remains challenging, with often no other solution than to resort to numerical simulations [28]. In this regard, two main types of DNA polymer models have been proposed. On one hand, explicit polymer models have allowed investigating supercoiling-induced phenomena on the basis of numerical simulations of single DNA chains [29, 30, 31, 32]. On the other hand, phenomenological models [33, 34] have led to bona fide analytical solutions of thermodynamic properties. These have been particularly useful to infer microscopic parameters of alternative DNA forms [25, 34], whose knowledge is crucial to both address the exact nature of these forms and parametrize explicit polymer models.

Explicit polymer models can be further divided into two types. Depending on the addressed question, structural details can indeed be coarse-grained at the scale of a single nucleotide [32, 35] or at a larger scale of a few tens base pairs [29, 30, 31, 36]. Single nucleotide resolution models have thus been useful to investigate properties of small (i.e. a few hundreds base pairs) DNA molecules as provided by cryo-electron microscopy [28, 37] and cyclization data [20], with the possibility to investigate sequence effects in detail (see e.g. [38] and references therein). Tens base pairs resolution models have instead been useful to address both mechanical and conformational properties of long (i.e. a few kb or tens kb) B-DNA molecules, taking advantage of much less time-consuming simulations and of the possibility to use parsimonious sets of parameters that can capture bending, torsional and self-avoidance properties of molecules [15, 20]. Along this line, the self-avoiding rod-like chain (sRLC) model [29] has been paramount to analyze folding properties of positively supercoiled molecules (see below for details of the model). These include molecular extensions [39], torques [40, 41, 42] and conformation details of superstructures [29, 43, 42] as measured by magnetic tweezers and cryo-electron microscopy, as well as dynamical properties (see [44] and references therein) and sequence-dependent phenomena [45, 46].

The sRLC model has thus provided a solid framework to infer physical parameters of B-DNA [15, 31]. To the best of our knowledge, it has however not yet been adapted to the co-existence of multiple DNA forms, precluding its use to investigate functional properties of bacterial genomes in vivo and leaving open several important biophysical questions. These include the balance between the global and local relaxations of supercoiling constraints, the exact nature of alternative forms (see e.g. [27] for a recent discussion about the experimental identification and characterization of denatured DNA (D-DNA) with respect to L-DNA) and, related to this, the inference of associated mechanical parameters. Namely, while some of the parameters, such as the free energy formation of left-handed DNA forms and of their associated junctions with B-DNA, have been estimated early on in bulk studies [22, 23, 47, 48, 49, 50, 51], the estimation of, e.g., the torsional and bending moduli of D-DNA and left-handed forms, has thus far relied exclusively on phenomenological models in the context of single molecule studies [33, 25, 26, 52, 34]. These models, which are based on bona fide free energy landscape descriptions of the co-existence of multiple DNA forms, have however never been quantitatively assessed, which may explain the existence of strong discrepancies in the prediction of some of the parameters [25, 34].

In this work, we aim at filling these gaps by extending the sRLC model so that it includes the possibility to have multiple DNA forms. To this end, we follow an approach initially proposed to tackle the problem of the large flexibility of small DNA molecules [53]. Namely, as proposed in [54, 55, 56], the softening of sharply bent DNA can be rationalized by considering a structural degree of freedom allowing the appearance of alternative forms more flexible than B-DNA, with a free energy cost reflecting the destabilization of the double helix. More elaborated models including torsional energies, but still neglecting DNA self-avoidance, have then been proposed to investigate the statistics of DNA denaturation [57] and the behavior of stretched DNA molecules [58] under torsional constraints. Models based on torsional energies alone have also been studied in the context of the statistics of DNA denaturation [59, 60] as well as to model properties of specific DNA sequences [26].

Following these studies, here we consider a version of the sRLC model, hereafter referred to as the 2sRLC model, which includes an additional structural degree of freedom. Using numerical simulations, we first demonstrate the ability of the 2sRLC model to reproduce magnetic tweezers experiments of long molecules (up to $21\text{\,}\mathrm{k}\mathrm{b}$) for a wide range of stretching forces and supercoiling levels typical of bacterial genomes in vivo. We next tackle the inference problem of parameters associated with alternative DNA forms. To this end, we circumvent the use of time-consuming simulations by deriving semi-analytical solutions of the thermodynamics when the writhe is negligible and we show that these can be used for any value of supercoiling level and stretching force where plectonemic superstructures are absent. As an application, we discuss the validity of the commonly used phenomenological model proposed by John Marko [33] and highlight a previously overlooked underdetermination problem related to the question of parameter inference, which is immanent to the high force analysis of single molecule extension curves. We provide, in this context, our best estimates of parameters associated with D-DNA and L-DNA.

2 Theory

2.1 The discrete sRLC model

In the sRLC framework, a double-stranded DNA molecule, usually in the form of B-DNA, is modeled as a continuous self-avoiding chain characterized by five fundamental parameters: i) a bending modulus ($K$) defining an associated persistence length ($\ell_{p}=K/k_{B}T$), ii) a torsional modulus ($C$) measured in units of twist persistence length, iii) a winding angle at rest ($\psi$) of the implicitly embedded helix, iv) the distance ($a$) between any two consecutive base pairs, thus defining the number of base pairs per $\ell_{p}$, and v) an effective radius ($r_{e}$) reflecting the hard-core simplification of the electrostatic repulsion of the DNA backbone [61]. For B-DNA in physiological conditions ($[\text{NaCl}]\sim$100\text{,}\mathrm{m}\mathrm{M}), $\ell_{p}$ typically lies in $[40-$60\text{\,}\mathrm{n}\mathrm{m}$]$ [[62](#bib.bib62), [63](#bib.bib63)] (in the following, we use $50\text{\,}\mathrm{n}\mathrm{m}$), $C$ is on the order of $100\text{\,}\mathrm{n}\mathrm{m}$ [[64](#bib.bib64), [24](#bib.bib24)], while $a=$0.34\text{\,}\mathrm{n}\mathrm{m}, $\psi=$0.6\text{,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{b}\mathrm{p} and $r_{e}\approx$2\text{\,}\mathrm{n}\mathrm{m} [61].

To simulate the folding of a DNA molecule, a discrete version of the model is used, for which the continuous formulation is a good approximation if the discretization is sufficiently fine (see hereafter for further details). Specifically, the chain is divided into identical impenetrable cylinders and an additional parameter, $n$, specifies the number of base pairs per cylinder used in the simulation such that the length of every cylinder is equal to $na$ (Fig. 1A) – in this work, we use $n\leq 10$ such that one cylinder corresponds at most to one B-DNA helix. The energy of any conformation $\mathcal{C}$ of a self-avoiding linear chain in the presence of a stretching force, $f$, then reads:

[TABLE]

where $z$ is the extension of the chain along the axis of the force, while $E_{B}(\mathcal{C})=\frac{k_{B}T}{2}\sum_{i=1}^{N}\frac{\ell_{p}}{na}\theta_{i}^{2}$ and $E_{T}(\mathcal{C})=\frac{k_{B}T}{2}\sum_{i=1}^{N}\frac{nC}{a}(\phi_{i}-\psi)^{2}$ respectively stand for the bending and torsional energies. In these equalities, the index $i$ indicates a site around which two consecutive cylinders are articulated, while $\theta_{i}$ is the bending angle associated with the $n$ base pairs of the site and $\phi_{i}$ is the winding angle per base pair (Fig. 1A); this means that $n\phi_{i}$ is the torsional counterpart of $\theta_{i}$ for the site $i$. Note that at the end of the molecule, two additional sites (at $i=0$ and $i=N+1$) set the boundary conditions, which here correspond to a situation where the orientation of the external cylinders is kept fixed along the force axis. Altogether, a linear chain made of $N_{\text{bp}}$ base pairs thus contains $N+1$ sites and $N$ segments ($N=N_{\text{bp}}/n=L/na$ where $L$ is the contour length of the molecule). A circular chain (e.g. a plasmid) would contain $N$ sites and $N$ segments, the first site joining the last and first segments together.

Importantly, just as $\theta_{i}$, $\phi_{i}$ can be explicitly written as a function of the local frames associated with the surrounding cylinders [65] (Fig. 1A). Using $\text{Tw}\equiv(2\pi)^{-1}\sum_{i=1}^{N}n\phi_{i}$, this ensures that the linking number, $\text{Lk}=\text{Tw}+\text{Wr}$, remains strictly constant for any deformation of a circular chain, provided the chain never crosses itself. Lk then reflects the topological status of the DNA molecule, for instance a supercoiling constraint if it is different from the corresponding value at rest, $\text{Lk}_{0}$. For a linear chain, conservation of the linking number further requires the ends of the chain to be attached to two impenetrable walls [39, 66]. These respectively play the role of the fixed surface and of the magnetic bead in single-molecule experiments.

2.2 The discrete 2sRLC model: including multiple DNA forms

In order to include varying DNA forms at every site $i$ of the sRLC, we follow the lines of denaturation studies [59, 54, 55, 56, 57, 58] and associate a state variable $s$ with every site ($s=\text{B, D, L or Z}$ if one considers for instance B-DNA, D-DNA, L-DNA and Z-DNA) specifying the associated mechanical parameters: $K_{s}$ (persistence length $\ell_{s}$), $C_{s}$, $a_{s}$ and $\psi_{s}$ – for simplicity, here we consider a form-independent electrostatic radius $r_{e}$. In this context, the length of cylinders may vary as a result of both a varying value for $a_{s}$ and a fixed number of base pairs per cylinder ($n$); here, for symmetry reasons, we consider the length of a cylinder surrounded by two sites with forms $s$ and $s^{\prime}$ to be equal to $n(a_{s}+a_{s^{\prime}})/2$ (Fig. 1B).

In addition to form-dependent mechanical parameters, alternative DNA forms come along with free energy formation costs on the order of one $k_{B}T$ per bp [67, 68], reflecting the internal deformation of base pairing and base pair stacking. Here, we denote this cost by $\gamma_{s}^{0}$ and set $\gamma_{B}^{0}=0$ (reference form). Finally, a domain wall penalty, $J$, must be considered between any two sites having different states. This penalty accounts for the internal free energy cost to transit from one DNA form to another one [69, 57, 52] (Fig. 1B) – in effect, it constrains alternative forms to gather in as few domains as possible. This wall penalty is in principle associated with two consecutive base pairs and therefore includes an entropic contribution on the order of $\ln(n)$ coming from the $n$ possible choices of these base pairs. For simplicity, here we consider this contribution as already included in the values of $J$.

Altogether, the energy of a self-avoiding 2sRLC conformation $\mathcal{C}$ having an arbitrary sequence $\{s_{i}\}_{i=1..N}$ of DNA forms reads:

[TABLE]

where $E_{\gamma}(\mathcal{C})=n\sum_{i=1}^{N}\gamma_{s_{i}}^{0}$ is the total free energy formation costs coming from alternative forms, $E_{B}(\mathcal{C})=\frac{k_{B}T}{2}\sum_{i=1}^{N}\frac{\ell_{s_{i}}}{na_{s_{i}}}\theta_{i}^{2}$ is the bending energy, $E_{T}(\mathcal{C})=\frac{k_{B}T}{2}\sum_{i=1}^{N}\frac{nC_{s_{i}}}{a_{s_{i}}}(\phi_{i}-\psi_{s_{i}})^{2}$ is the torsional energy and $E_{J}(\mathcal{C})=J\sum_{i=1}^{N-1}(1-\delta_{s_{i},s_{i+1}})$ is the total wall penalty, with $\delta_{s_{i},s_{i+1}}=1$ if $s_{i}=s_{i+1}$, [math] otherwise.

3 Results and discussion

3.1 Effective 2sRLC model

Solving the thermodynamics of the 2sRLC model is challenging because of the self-avoidance constraints. Yet, following upon previous studies on single DNA forms [70] and on the statistics of DNA denaturation [59, 60], it is possible to integrate out the torsional degrees of freedom (the $\phi_{i}$’s) under the constraint of a fixed Lk. This leads to an equivalent effective model that has the benefit to offer much better simulation performances [39, 42] and further analytic treatment (see below). In the case of the 2sRLC model, the resulting conformational energy can be decomposed as in Eq. 2, but with an effective free energy of formation and an effective torsional energy that respectively read (see Supplementary methods inSupplementary material for details):

[TABLE]

where, again, $N_{\text{bp}}=nN$ is the length of the molecule in base pairs, \psi(\mathcal{C})\equiv\sum_{i=1}^{N}\psi_{s_{i}}\big{/}N and {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}(\mathcal{C})\equiv\sum_{i=1}^{N}\frac{a_{s_{i}}}{C_{s_{i}}}\big{/}N are conformation-dependent average quantities, namely, the average twist angle per bp at rest and the average torsional susceptibility per bp.

Eq. 4 thus generalizes to the case of any number of different DNA forms the possibility, for topologically constrained DNA molecules, to obtain an equivalent effective model with a torsional energy quadratic in Tw. Note, nevertheless, that here the associated parameters depend on the composition of the molecule in these DNA forms. In the particular case of a single (B-DNA) form, one for instance recovers $\psi(\mathcal{C})=\psi_{B}$ and {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}(\mathcal{C})=a_{B}/C_{B} for all conformations, such that $E_{T}^{eff}(\mathcal{C})=\frac{k_{B}TC_{B}N_{\text{bp}}}{2a_{B}}\left(2\pi\frac{\text{Tw}(\mathcal{C})}{N_{\text{bp}}}-\psi_{B}\right)^{2}$ as demonstrated in [70].

3.2 Capturing the phenomenology of negatively supercoiled DNA

To demonstrate the ability of the 2sRLC model to capture the phenomenology of negatively supercoiled DNA, we performed Monte-Carlo simulations of long chains using the above effective energies. To this end, we have adapted the standard Monte-Carlo procedure used to simulate the equilibrium folding of the sRLC model [31] in order to include the possibility for mechanical parameters to change at any site of the chain, at any iteration of the simulation (see [66] for details). In practice, it involves a new transition type so that for any site $i$ independently of the other sites, the value of the variable $s_{i}$ may change under the constraint of detailed balance. Corresponding variations of cylinder lengths (Fig. 1B) are thus managed as other elementary Monte-Carlo moves. In particular, if they lead to a collision, the transition is rejected. More generally, length variation of a cylinder always requires an adjustment of the positions and orientations of the surrounding cylinders to preserve the continuity of the chain. As explained in [66], this is realised by performing specific random rotations that are constrained by the geometry of the chain.

Using this Monte-Carlo method, we could quantitatively reproduce extensions of several kilo base pairs long molecules as obtained in magnetic tweezers experiments for a wide range of forces and supercoiling densities $\sigma=(\text{Lk}-\text{Lk}_{0})/\text{Lk}_{0}$. As an example, we report in Fig. 2 both our simulation results and the experimental results of Cees Dekker’s lab [27] for a 21 kb molecule, for $\sigma\in[-0.05,0.05]$, $f=$0.5\text{,}\mathrm{p}\mathrm{N} and $f=$4.5\text{\,}\mathrm{p}\mathrm{N}, providing as a bonus a prediction of the equilibrium fraction of alternative forms present along the molecule (top panel). DNA denaturation being expected to be the most likely alternative form at these forces and supercoiling densities [26, 27], here we considered the situation where only B-DNA and D-DNA could form, with the following parameters leading to good agreement with experimental curves: $a_{D}=$0.45\text{,}\mathrm{n}\mathrm{m}, $\ell_{D}=$15\text{\,}\mathrm{n}\mathrm{m}, $C_{D}=$10\text{,}\mathrm{n}\mathrm{m} and $\gamma_{D}=2k_{B}T$ (see below for discussion of these values) – we intentionally used a $\ell_{D}$ much larger than the $4\text{\,}\mathrm{n}\mathrm{m}$ diameter of the chain because smaller values would lead, in any case, to larger effective persistence lengths [[71](#bib.bib71)]. We used $J=10k_{B}T$, in accord with previous analyses [[69](#bib.bib69), [26](#bib.bib26), [52](#bib.bib52)], and considered a discretization level of $n=10$, leading to $\simeq 3.3$ cylinders per $\ell_{D}$ and 15 cylinders per $\ell_{B}(=$50\text{\,}\mathrm{n}\mathrm{m}$)$, a value that both offers reasonable computational times and ensures properties of B-DNA superstructuring to be insensitive to discretization procedures [[42](#bib.bib42), [66](#bib.bib66)]. Finally, the quantitative reproduction of the buckling transition at $\sigma>0$ (right part of the curves in Fig. [2](#S3.F2)) required us to use $C=$80\text{\,}\mathrm{n}\mathrm{m} at $f=$4.5\text{,}\mathrm{p}\mathrm{N} and $C=$100\text{\,}\mathrm{n}\mathrm{m} at $f=$0.5\text{,}\mathrm{p}\mathrm{N}$$ (Fig. S1). This is in accord with the observed systematic discrepancies occurring at low forces between sRLC and experimental extensions [72], a consequence of a coupling between torsional and bending deformations stemming from the asymmetry of the DNA helix [73].

Our simulations corroborate well-known results of magnetic tweezers experiments, with notably the presence of two regimes (Fig. 2): i) a low force regime ($f\lesssim$0.5\text{,}\mathrm{p}\mathrm{N}) that is symmetric under the transformation $\sigma\to-\sigma$ and where torsional stress is partially relaxed through the formation of plectonemes, some of them being branched, and ii) a $\sigma$-asymmetric high force regime ($f\gtrsim$2\text{\,}\mathrm{p}\mathrm{N}) where negative torsional stress is partially relaxed by forming denaturation bubbles. For intermediate forces at $\sigma<0$, both conformations with denaturation bubbles and plectonemic conformations become locally stable, being separated by high free energy barriers (Fig. S2). This is in accord with the experimental findings, in this regime, of large equilibrium fluctuations of the extension  [74, 34], which have been shown, in the case of a 2.4 kb long molecule, to reflect a two-state behavior between plectonemic conformations and conformations with D-DNA bubbles [34]. In our simulations of comparable 2.4 kb long molecules, this manifests through an hysteresis cycle in the force-extension diagram (Fig. S2), with the drawback that specific numerical methods, such as metadynamics [75], must be developed in this regime to get access to thermodynamic properties, including the height of the free energy barriers.

Our simulations also corroborate the existence of conformations where denaturation bubbles locate at the apex of plectonemes (Fig. S3). These have been predicted to occur using a polymer model of DNA coarse-grained at the nucleotide level [76]. Here, we observe that these tip-bubble conformations, as coined in [76], are all the more favored that the system is close to the ”spinodal” point associated with the plectonemic state – the force at which the free energy barrier separating pure plectonemic states from states with denaturation bubbles becomes on the order of $k_{B}T$.

Most importantly, while trying to reproduce previously published experimental extensions, we have found that different sets of parameters could lead to similar results. For instance, the D-DNA persistence length used in Fig. 2 is typically 5 times larger than that previously estimated  [33, 25], while previous works have led to values of the torsional modulus that may differ by more than 20-fold [25, 34]. Overall, this raises the question of the extent to which it is possible to infer mechanical parameters of alternative forms using extension curves alone as obtained in magnetic tweezers experiments, what we discuss now in detail.

3.2.1 Zero writhe approximation in the high force regime

One could in principle use the above numerical simulations to search for sets of parameters that best reproduce extension curves. However, independently of the efficiency of simulations, an exhaustive exploration of the 5 parameters ($\ell_{D}$, $C_{D}$, $a_{D}$, $\psi_{D}$ and $\gamma_{D}^{0}$) and the domain wall penalty $J$ is extremely challenging. Moreover, the high free energy barriers separating plectonemic conformations from conformations with D-DNA bubbles precludes obtaining the relative fraction of each state at equilibrium using standard Monte-Carlo methods. Curves at intermediate forces between $0.5\text{\,}\mathrm{p}\mathrm{N}$ and $2\text{\,}\mathrm{p}\mathrm{N}$ cannot thus be currently exploited to further constrain parameters as previously done in the context of phenomenological models [34].

To partially circumvent these problems, we focus on the regime of high forces, for which the thermodynamics can be solved in the approximation of both a negligible writhe, i.e. when $\text{Lk}=\text{Tw}$, and negligible self-avoidance properties (see below for the range of forces where this approximation holds). In this approximation, it is indeed possible to obtain an explicit formula for the free energy profile of the system, $F_{f,\sigma}(\{\lambda_{i}\})$, as a function of the fraction of each of the possible alternative forms at a given stretching force $f$ and supercoiling density $\sigma$. In the following, for simplicity we discuss the case of two forms, with in mind B-DNA plus a certain fraction $\lambda$ of an alternative form, hereafter denoted as X-DNA – the generalization to more than two forms is straightforward (see Supplementary methods).

Calling $\mathcal{J}(\lambda)$ the free energy associated with the multiple possibilities to distribute the $\lambda N$ sites of the alternative form into distinct domains (see SI for its derivation), $g_{s,f}(L_{s})$ the free energy of the WLC at force $f$ associated with the form $s$ with contour length $L_{s}$, and using the notations {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{s}=a_{s}/C_{s} and {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}/\psi(\lambda)=(1-\lambda){\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{B}/\psi_{B}+\lambda{\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{s}/\psi_{s}, we obtain (Supplementary methods):

[TABLE]

Using this free energy profile, equilibrium properties of the 2sRLC model in the zero writhe approximation can be computed by performing numerical integration over $\lambda$ of Boltzmann-weighted observables (see Eq. 10 in Supplementary methods and subsequent explanations) with, here, a particular interest in equilibrium values of the extension, $z^{*}=-\langle\partial_{f}F_{f,\sigma}(\lambda)\rangle_{\lambda}$, the torque $\Gamma^{*}=(2\pi\text{Tw}_{0})^{-1}\langle\partial_{\sigma}F_{f,\sigma}(\lambda)\rangle_{\lambda}$ (see Supplementary methods for derivation and further details), and the fraction of alternative form, $\lambda^{*}=\langle\lambda\rangle_{\lambda}$, where \langle\bullet\rangle_{\lambda}\equiv\int d\lambda\bullet\exp\left[-F_{f,\sigma}(\lambda)/k_{B}T\right]\big{/}\int d\lambda\exp\left[-F_{f,\sigma}(\lambda)/k_{B}T\right]. We also compute the equilibrium number of domains, $X^{*}$, using a free energy surface in the $(\lambda,X)$ space (Eq. 6 in Supplementary methods).

3.2.2 Validity domain of the zero writhe approximation

To assess the validity domain of the zero-writhe approximation, we compare the above Boltzmann-weighted equilibrium values to those obtained in Monte-Carlo simulations of the corresponding 2sRLC model. As a case study, we investigate a 1 kb long molecule with an alternative form characterized by $\ell_{X}=$15\text{,}\mathrm{n}\mathrm{m}, $a_{X}=$0.54\text{\,}\mathrm{n}\mathrm{m}, $\psi_{X}=0$, $C_{X}=$10\text{,}\mathrm{n}\mathrm{m} and $\gamma_{X}=1.6kT$, having in mind a form close to D-DNA; just as in the analysis of Fig. [2](#S3.F2), the relatively high value of $\ell_{X}$ compared to previous estimations for D-DNA ($3-$4\text{\,}\mathrm{n}\mathrm{m} [33, 25]) was chosen for practical purposes, that is, to avoid potential artifacts that may occur in the simulations when the persistence length is smaller or on the order of $r_{e}$ [66] and to avoid considering (longer) effective persistence lengths in our analytical approximations [71]. We also used a fine discretization level of $n=4$ to get a large number of cylinders per $\ell_{D}$ ($\simeq 7$) so that to prevent any discrepancy coming from the discretization procedure of the rod-like chain.

Regarding $\lambda^{*}$ and the extension, we observe an excellent agreement between simulation results and the zero-writhe approximation for any force above $2\text{\,}\mathrm{p}\mathrm{N}$ (left panels in Fig. 3), with a relative error lower than 5% (Fig. S4). For the number of domains the error is on the order of 20% (a higher error is expected as a consequence of the discrete nature of $X$), and up to 30% for the torque (side effect of the error in $X^{*}$). Most remarkably, this agreement is observed for a wide range of supercoiling densities, from the transition point at physiological $\sigma\simeq-0.025$ associated with the onset of the first denaturation bubble (see inset in upper right panel of Fig. 3) to large non-physiological negative values of $\sigma=-0.5$, and for forces down to the onset of plectonemic superstructures. The strong deviation observed at $1\text{\,}\mathrm{p}\mathrm{N}$ is indeed the result of a buckling transition that occurs in the simulations, which is not captured by the zero writhe approximation. In other words, the zero writhe approximation provides a good description of the system whenever plectonemic superstructures are absent, paving the way for an efficient method for the quantitative inference of physical properties of alternative DNA forms in the high force regime (see below). In this regard, we note that the inclusion of finite size corrections [77] for computing the WLC terms $g_{X,f}$ and $g_{B,f}$ leads to very similar results (Fig. S5), suggesting the possibility to use the classical infinite length approximation of the WLC [78], although the fraction of alternative form may be large.

We finally used the same simulations to assess the validity of the commonly used phenomenological approach proposed by John Marko a decade ago [33]. In a nutshell, Marko proposed a model with two separated phases (B-DNA and X-DNA), each of them being allocated with a certain supercoiling density and a free energy that is quadratic in the latter; the two supercoiling densities are coupled because of the conservation of the linking number. Then, supposing a linear relationship between the fraction $\lambda$ of X-DNA and the total supercoiling density ($\sigma$) and considering the torques to be equal in both phases, Marko used a free energy minimization procedure to compute the supercoiling densities specific to each phase, leading back to the computation of the free energy as a function of $\sigma$. Here, our analysis shows that this phenomenological approach provides as well a remarkable description of the 2sRLC equilibrium properties (Fig. 3), excluding the number of domains, $X$, which is not a variable of the model [33].

3.3 Inferring mechanical parameters of alternative forms from extension curves

Alternative DNA forms mostly occur in presence of B-DNA. Estimation of their parameters, which ultimately reflect statistical properties of the dynamics of strands, has thus been a matter of debate. Both bulk [67] and single-molecule [68] measurements have nevertheless converged to free energy formations of D-DNA on the order of one $k_{B}T$ per bp. Cyclization measurements of small circular DNA molecules have also revealed an independence of D-DNA strands, suggesting to use $\phi_{D}=0$ [79]. A zero-writhe phenomenological modeling of these experiments have then led to $a_{D}=$0.54\text{,}\mathrm{n}\mathrm{m}, $\ell_{D}\simeq 2\text{ to }3~{}\text{nm}$ and $C_{D}=$1\text{\,}\mathrm{n}\mathrm{m} [79]. In contrast, investigation of the co-existence of D-DNA with B-DNA at forces $\lesssim$1\text{,}\mathrm{p}\mathrm{N} and $|\sigma|\lesssim 0.06$, using a phenomenological model of the coexistence of twisted B-DNA, plectonemic B-DNA, and D-DNA, has led to $C_{D}=$28\text{\,}\mathrm{n}\mathrm{m} [34]. Regarding left-handed DNA forms, static and dynamic laser light scattering experiments of Z-DNA solutions led to $\ell_{Z}\simeq$200\text{,}\mathrm{n}\mathrm{m} [[80](#bib.bib80)], whereas a best fitting procedure using Marko’s phenomenological model revealed a flexible L-DNA form, with $a_{L}=$0.48\text{\,}\mathrm{n}\mathrm{m} and $\ell_{L}=$3\text{,}\mathrm{n}\mathrm{m} [[25](#bib.bib25)], and a torsional modulus $C_{L}$ lying between $10$ and $20\text{\,}\mathrm{n}\mathrm{m}$. Similar values of $C_{L}$ were found by investigating the collapse of “pure” L-DNA molecules  [[26](#bib.bib26)]. Finally, location of the buckling transition of such ”pure” L-DNA molecules at high forces led to $\psi_{L}=$-0.5\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{b}\mathrm{p} [26, 27].

Here, we aim at using the predictive power of our zero writhe approximation of the 2sRLC model (Fig. 3) in order to infer some of the parameters of alternative forms. Specifically, we discuss the possibility, as proposed in [25], to estimate parameters of D-DNA and L-DNA using $\sigma$-extension curves, alone, in the high force regime. To this end, we first test our capacity to recover original parameters using simulated data and, then, discuss estimations from real data.

3.3.1 Simulated data: assessing parameter inference in an underdetermination context

To test our capacity to recover original parameters (hereafter indicated by a hat, e.g., $\hat{a}_{X}$) from simulated $\sigma$-extension data, we use as a generative model the 2sRLC model under the zero writhe approximation (Eq. 5) and question the possibility to recover original parameters using the same generative model. Regarding the inference method, we follow the procedure of [25], which consists in determining optimal parameters such that extensions of the model best fit experimental extensions over a wide range of forces and supercoiling densities (Supplementary methods).

First, compared to the other parameters, we find that $a_{X}$ and $\ell_{X}$ (or equivalently $K_{X}$) have a more systematic impact on $\sigma$-extension curves. More precisely, using a large random set of parameters, the various $\sigma$-extension curves appear to be sorted according to the values of both $a_{X}$ and $\ell_{X}$, whereas they appear to be randomized with respect to the values of the other parameters (Fig. 4A, Fig. S6) – we nevertheless note that best fits correspond to values of $\psi_{X}$ that are distant from $\psi_{B}$ (Fig. S6). The goodness of fit can be further quantified by computing the root-mean-square deviation (RMSD) of each tested curve with respect to the original curve. Reporting these RMSDs in the plane $(a_{X},\ell_{X})$ at a given force then reveals the existence of a well-defined continuous set of parameters with similar locally optimal solutions and going through $(\hat{a}_{X},\hat{\ell}_{X})$, hereafter called a crest of solutions (Fig. 4B) – note that these solutions are always suboptimal with respect to the original set of parameters. The overall shape of this crest reflects the existence of two types of solutions: 1) solutions corresponding to almost straight pieces of X-DNA (large $\ell_{X}$, the exact value being irrelevant); 2) solutions corresponding to more disordered domains (small $\ell_{X}$) with longer double helices (larger $a_{X}$). Most importantly, an analysis of solutions in the vicinity of $(\hat{a}_{X},\hat{\ell}_{X})$ corroborates the existence of a large spectrum of possible values for $J$, $\gamma_{X}$, $\psi_{X}$ and $C_{X}$ (Fig. S8), meaning that these parameters cannot be estimated using a single extension curve, i.e., a single force. Moreover, the presence of the crest implies the existence of some degeneracy that precludes the unambiguous determination of $\hat{a}_{X}$ and $\hat{\ell}_{X}$, too.

In principle, the value of $\hat{C}_{X}$ could be addressed using torque measurements [26], yet only if one has access to $\lambda(\sigma)$ (see Supplementary methods). For $\hat{a}_{X}$ and $\hat{\ell}_{X}$, one might expect that the underdetermination problem akin to a single force could be resolved by considering multiple forces. However, we find that, without fixing any other parameters, the crest gets blurrier as the force increases (blue points in Fig. 4C). To understand this intriguing result, we performed the same analyses but by setting one arbitrary parameter to its right value. We found that fixing either $J$ or $\gamma_{X}$ leads to the same underdetermination problem. In contrast, fixing either $\psi_{X}$ or $C_{X}$ makes a crest of locally suboptimal solutions re-appear. However, setting $\psi_{X}$ to the right value $\hat{\psi}_{X}=0$ results in a crest that is off the right solution $(\hat{a}_{X},\hat{\ell}_{X})$ (Fig. S9), a shift that disappears only by setting $C_{X}=\hat{C}_{X}$ (red points in Fig. 4C).

Overall, these results show that inference of $\hat{a}_{X}$ and $\hat{\ell}_{X}$ from extension curves requires to have some prior knowledge about the nature of the alternative form, more particularly about its torsional parameters ($\hat{\psi}_{X}$ and $\hat{C}_{X}$). This may be explained by the fact that the (effective) torsional energy is quadratic in $\sigma$, such that any parameter that affects this energy is expected to have a large impact on any curve plotted as a function of $\sigma$. When the values of $\hat{\psi}_{X}$ and $\hat{C}_{X}$ are known, a minimization over all forces then lead to a fair estimation of both $\hat{a}_{X}$ and $\hat{\ell}_{X}$ which is all the more accurate that the range of available forces is large (the crossing of the different crests of solutions becoming well-defined as shown in Fig. S7).

3.3.2 Experimental data: estimation of D-DNA and L-DNA parameters

To estimate $\hat{a}$ and $\hat{\ell}_{p}$ for both D-DNA and L-DNA, we use the data of Sheinin et al. [25] obtained for 2.2 kb long molecules, which consist of extension curves measured over a wide range of values of $\sigma$ (from -2.3 to +0.1) and $f$ (from 2.5 to $36\text{\,}\mathrm{p}\mathrm{N}$). In accord with the findings of [27], we further consider that L-DNA occurs for $\sigma\lesssim-1$ and, hence, estimate D-DNA and L-DNA parameters using an optimization procedure of the extension curves for $\sigma$ lying in $[-0.6,-0.1]$ and in $[-1.8,-1.3]$, respectively. As discussed above, we lift the underdetermination problem by setting, on one hand, $\psi_{D}=0$ and $\psi_{L}=$-0.5\text{,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{b}\mathrm{p} and we test, on the other hand, $C_{D}=$1\text{\,}\mathrm{n}\mathrm{m} as discussed in [25], $C_{D}=$28\text{,}\mathrm{n}\mathrm{m} as proposed in [[34](#bib.bib34)] and $C_{L}=$10\text{\,}\mathrm{n}\mathrm{m} in accord with previous estimations [25, 26].

Similar to the case of simulated data, optimizations over forces from $2.5$ to $12\text{\,}\mathrm{p}\mathrm{N}$ reveals the existence of small islands of optimal solutions, which we consider as our best estimations of parameters (Fig. 5). For D-DNA (Fig. 5A), $C_{D}=$1\text{,}\mathrm{n}\mathrm{m} leads to $a_{D}>$0.6\text{\,}\mathrm{n}\mathrm{m} and $\ell_{D}$ close to $1\text{\,}\mathrm{n}\mathrm{m}$, which are different, although being qualitatively similar, from previous estimations according to which $a_{D}=$0.54\text{,}\mathrm{n}\mathrm{m} and $\ell_{D}\in[2-$3\text{\,}\mathrm{n}\mathrm{m}$]$ (see [[25](#bib.bib25)] and references therein). When $C_{D}=$28\text{\,}\mathrm{n}\mathrm{m}, we find $a_{D}\in[0.45-$0.5\text{,}\mathrm{n}\mathrm{m}$]$ and $\ell_{D}\in[2-$5\text{,}\mathrm{n}\mathrm{m}$]$. Note, here, that we used the finest level of chain discretization ($n=1$) and we checked that our results remain identical for rougher coarse-graining up to $n=4$ (Fig. S10). The effect of $n$ is indeed expected to affect the exact value of the domain wall penalty $J$ (see above) and the entropic contribution to the free energy associated with the multiple possible positions of the denaturation bubbles (the term $\mathcal{J}(\lambda)$ in Eq. 5, see Supplementary methods). However, as discussed using simulated data, the former is expected to play only a marginal role in the determination of $\hat{a}_{D}$ and $\hat{\ell}_{D}$, while in the case of a fixed, small number of denaturation bubbles as here (typically one), the latter scales as $\ln(n)$.

For L-DNA (Fig. 5B), we find $a_{L}\in[0.5-$0.6\text{,}\mathrm{n}\mathrm{m}$]$, which is slightly larger than previous estimations [25], and $\ell_{L}\in[1-$2\text{,}\mathrm{n}\mathrm{m}$]$, which is slightly smaller. In all cases, the small value of $\ell_{L}$ corroborates that L-DNA is not a pure left-handed structural form but made of denatured DNA [25, 26].

To summarize, our results question the compatibility of previously used $a_{D}=$0.54\text{,}\mathrm{n}\mathrm{m}$$ with reported torsional parameters. They also show the importance of having precise knowledge of a certain number of parameters to be able to fully exploit magnetic tweezers experiments in the high force regime. In absence of such prior information (as it is usually the case for real molecules), several strategies can be contemplated. For instance, working with an optimization procedure that is constrained by both torque curves and extension curves should narrow down the set of parameters leading to good fits. Even in this case, yet, having access to the fraction of the alternative form seems mandatory (see Supplementary methods). Potential development along this line might come from the use of DNA binding proteins sensitive to alternative forms. Additional constraints on model parameters could also be imposed using cyclization data and conformation statistics from cryo-electron microscopy experiments.

3.3.3 Extension of the model by including sequence effects

In effect, our 2sRLC model bridges ”kinkable” worm-like chain models and self-avoiding rod like chain models, which have been developed to respectively address the high flexibility of sharply bent DNA [53] and the superstructuring properties of supercoiled B-DNA molecules [15, 20]. By doing so, it opens novel perspectives in the field of supercoiled DNA, not only to better rationalize magnetic tweezers experiments but also to predict behaviors of bacterial genomes in vivo. In this regard, including sequence effects, such as the tendency of AT-tracks and CpG-tracks to respectively favor DNA denaturation and left-handed forms [81], can be easily realized by using sequence dependent free energy formation costs of alternative forms (the $\gamma_{s_{i}}$ in Eq. 3). As an example, we simulated the folding of a 1 kb long molecule where one single site ($i_{AT}$) had a much lower $\gamma^{0}_{D}$ than the other sites, thus mimicking the presence of an AT-track. Just as previously found in a more detailed, base resolution polymer model of DNA [76], we find that for a force $f=$1\text{,}\mathrm{p}\mathrm{N}$$ and a supercoiling density $\sigma=-0.06$, the site $i_{AT}$ is systematically denatured and serves as both nucleation and anchor points for a plectoneme (Fig. S11).

4 Conclusion

More elaborate 2sRLC models can be contemplated by e.g. including sequence effects at the level of DNA bending [82, 38]. In particular, by properly distinguishing intrinsic curvature from stiffness [38], one can expect to quantitatively address the structuring properties of bacterial sequences, including the possible pinning of pure B-DNA plectonemes [83]. In all cases, our ability to develop truly predictive models of supercoiled DNA depends on our capacity to precisely infer mechanical parameters of alternative forms, which itself relies on the derivation of analytical expressions of thermodynamic quantities as a function of system parameters, at least for some range of forces and supercoiling densities [25, 26, 34]. Here, we have derived simple analytical expressions that accurately describe the free energy landscape of the 2sRLC model in the approximation of zero writhe, providing a powerful alternative to phenomenological approaches. It would then be interesting to compare our analytical results to other exact analytical results obtained in models devoid of self-avoidance [58]. It would also be interesting to further expand the free energy landscape in power series of the writhe. This would allow including additional constraints coming from data at lower forces where alternative forms co-exist with superstructures [34].

This latter possibility appears to be crucial for unambiguously estimating parameters of alternative DNA forms. We have indeed shown that the proper estimation of a specific parameter may critically depend on the knowledge of other parameters or, equivalently, on their independent estimation – see e.g. the impact of the torsional modulus on the estimation of both the distance between consecutive base pairs and the persistence length in Fig. 5B. A systematic parallel analysis of extension curves together with torque curves (which is particularly adapted to highlight torsional properties) seems therefore necessary. In addition, sequences of interest, like GC or AAT repeats as implemented in  [26, 52], could be systematically used to channel the formation of alternative forms, which would have two main advantages. On one hand, this would allow ”controlling” the fraction of the alternative form ($\lambda$ in Eq. 5), thus lifting an important unknown of the problem. On the other hand, this would sharpen signals because of stronger cooperative transitions (see e.g. [52]) and, hence, should constrain further model predictions – just as for the co-existence of plectonemic conformations and strecthed conformations with alternative forms (Figs. 2 and S2), this may also require specific numerical methods to efficiently sample the different states at equilibrium. For more general sequences, DNA binding proteins sensitive to alternative forms, as used e.g. in [27], might offer solutions to compute the fraction $\lambda$.

Acknowledgements

We thank Ruggero Cortini and Marc Joyeux for their early input in this project, Bahram Houchmandzadeh for useful suggestions, Daniel Jost for critical reading of the manuscript and Olivier Rivoire for helpful comments. We also thank Cees Dekker and John van Noort for providing us with their experimental data. I.J. is supported by an ATIP-Avenir grant (Centre National de la Recherche Scientifique).

Appendix A. Supplementary Material

The supplementary material is a single pdf file consisting of two sections: I) the ”Supplementary methods” section provides details related to the analytics and simulation of the 2sRLC model and II) the ”Supplementary figures” section contains 11 figures.

Author Contributions

I.J and T.L. performed research; T.L. implemented simulation tools; I.J. and T.L. analyzed data; I.J. wrote the article.