Thibaut Lepage

Ivan Junier

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

Keywords: DNA polymer modeling, DNA supercoiling, Multiscale polymer models, DNA structural transitions

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.

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 [START_REF] Liu | Supercoiling of the DNA template during transcription[END_REF][START_REF] Postow | Topological challenges to DNA replication: conformations at the fork[END_REF]. 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 [START_REF] Junier | Universal and idiosyncratic characteristic lengths in bacterial genomes[END_REF], this negative supercoiling is further partitioned into approximately 10 kb long topologically independent domains [START_REF] Postow | Topological domain structure of the Escherichia coli chromosome[END_REF][START_REF] Deng | Organization of supercoil domains and their reorganization by transcription[END_REF].

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 [START_REF] Strick | Stretching of macromolecules and proteins[END_REF]. Thus, for topologically constrained DNA as in the case of circular molecules (e.g. plasmids), of topologically constrained linear domains [START_REF] Liu | Supercoiling of the DNA template during transcription[END_REF], or when DNA molecules are manipulated by magnetic tweezers [START_REF] Bryant | Recent developments in singlemolecule DNA mechanics[END_REF][START_REF] Lionnet | Single-molecule studies using magnetic traps[END_REF][START_REF] Kriegel | Probing the mechanical properties, conformational changes, and interactions of nucleic acids with magnetic tweezers[END_REF], 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 [START_REF] Du | The genome-wide distribution of non-B DNA motifs is shaped by operon structure and suggests the transcriptional importance of non-B DNA structures in Escherichia coli[END_REF][START_REF] Vologodskii | Biophysics of DNA[END_REF], including denaturation bubbles [START_REF] Benham | Torsional stress and local denaturation in supercoiled DNA[END_REF], cruciforms [START_REF] Benham | Stable cruciform formation at inverted repeat sequences in supercoiled DNA[END_REF], left-handed Z-DNA [START_REF] Rich | The chemistry and biology of lefthanded Z-DNA[END_REF] and lefthanded L-DNA [START_REF] Bryant | Structural transitions and elasticity from torque measurements on DNA[END_REF][START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF].

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 [START_REF] Irobalieva | Structural diversity of supercoiled DNA[END_REF]. 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 [START_REF] Vologodskii | Conformational and thermodynamic properties of supercoiled DNA[END_REF][START_REF] Klenin | A Brownian dynamics program for the simulation of linear and circular DNA and other wormlike chain polyelectrolytes[END_REF][START_REF] Vologodskii | Simulation of equilibrium and dynamic properties of large DNA molecules, Computational Studies of DNA and RNA[END_REF][START_REF] Ouldridge | Structural, mechanical, and thermodynamic properties of a coarse-grained DNA model[END_REF]. On the other hand, phenomenological models [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF] have led to bona fide analytical solutions of thermodynamic properties. These have been particularly useful to infer microscopic parameters of alternative DNA forms [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF], 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 [START_REF] Ouldridge | Structural, mechanical, and thermodynamic properties of a coarse-grained DNA model[END_REF][START_REF] Manghi | Physics of base-pairing dynamics in DNA[END_REF] or at a larger scale of a few tens base pairs [START_REF] Vologodskii | Conformational and thermodynamic properties of supercoiled DNA[END_REF][START_REF] Klenin | A Brownian dynamics program for the simulation of linear and circular DNA and other wormlike chain polyelectrolytes[END_REF][START_REF] Vologodskii | Simulation of equilibrium and dynamic properties of large DNA molecules, Computational Studies of DNA and RNA[END_REF][START_REF] Krajina | Large-Scale Conformational Transitions in Supercoiled DNA Revealed by Coarse-Grained Simulation[END_REF]. 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 [START_REF] Irobalieva | Structural diversity of supercoiled DNA[END_REF][START_REF] Wang | Influence of DNA sequence on the structure of minicircles under torsional stress[END_REF] and cyclization data [START_REF] Vologodskii | Biophysics of DNA[END_REF], with the possibility to investigate sequence effects in detail (see e.g. [START_REF] Mitchell | Maddocks, Sequence-Dependent Persistence Lengths of DNA[END_REF] 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 [START_REF] Strick | Stretching of macromolecules and proteins[END_REF][START_REF] Vologodskii | Biophysics of DNA[END_REF]. Along this line, the self-avoiding rod-like chain (sRLC) model [START_REF] Vologodskii | Conformational and thermodynamic properties of supercoiled DNA[END_REF] has been paramount to analyze folding properties of positively supercoiled molecules (see below for details of the model). These include molecular extensions [START_REF] Vologodskii | Extension of torsionally stressed DNA by external force[END_REF], torques [START_REF] Lipfert | Magnetic torque tweezers: measuring torsional stiffness in DNA and RecA-DNA filaments[END_REF][START_REF] Schöpflin | Probing the Elasticity of DNA on Short Length Scales by Modeling Supercoiling under Tension[END_REF][START_REF] Lepage | Thermodynamics of Long Supercoiled Molecules: Insights from Highly Efficient Monte Carlo Simulations[END_REF] and conformation details of superstructures [START_REF] Vologodskii | Conformational and thermodynamic properties of supercoiled DNA[END_REF][START_REF] Bednar | The twist, writhe and overall shape of supercoiled DNA change during counterion-induced transition from a loosely to a tightly interwound superhelix. Possible implications for DNA structure in vivo[END_REF][START_REF] Lepage | Thermodynamics of Long Supercoiled Molecules: Insights from Highly Efficient Monte Carlo Simulations[END_REF] as measured by magnetic tweezers and cryo-electron microscopy, as well as dynamical properties (see [START_REF] Ivenso | Simulation of DNA Supercoil Relaxation[END_REF] and references therein) and sequence-dependent phenomena [START_REF] Klenin | Modulation of intramolecular interactions in superhelical DNA by curved sequences: a Monte Carlo simulation study[END_REF][START_REF] Chirico | Brownian dynamics simulations of supercoiled DNA with bent sequences[END_REF].

The sRLC model has thus provided a solid framework to infer physical parameters of B-DNA [START_REF] Strick | Stretching of macromolecules and proteins[END_REF][START_REF] Vologodskii | Simulation of equilibrium and dynamic properties of large DNA molecules, Computational Studies of DNA and RNA[END_REF]. 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. [START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF] 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 [START_REF] Benham | Stable cruciform formation at inverted repeat sequences in supercoiled DNA[END_REF][START_REF] Rich | The chemistry and biology of lefthanded Z-DNA[END_REF][START_REF] Singleton | Left-handed Z-DNA is induced by supercoiling in physiological ionic conditions[END_REF][START_REF] Peck | Energetics of B-to-Z transition in DNA[END_REF][START_REF] Haniford | Facile transition of poly[d(TG) x d(CA)] into a left-handed helix in physiological conditions[END_REF][START_REF] Vologodskii | Left-handed z form in superhelical dna: A theoretical study[END_REF][START_REF] Bauer | The free energy, enthalpy and entropy of native and of partially denatured closed circular DNA[END_REF], the estimation of, e.g., the torsional and bending moduli of D-DNA and lefthanded forms, has thus far relied exclusively on phenomenological models in the context of single molecule studies [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF][START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Oberstrass | Torque Spectroscopy of DNA: Base-Pair Stability, Boundary Effects, Backbending, and Breathing Dynamics[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. 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 [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF].

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 [START_REF] Cong | Recent progress on the mechanics of sharply bent DNA[END_REF]. Namely, as proposed in [START_REF] Yan | Localized single-stranded bubble mechanism for cyclization of short double helix DNA[END_REF][START_REF] Wiggins | Exact theory of kinkable elastic polymers[END_REF][START_REF] Ranjith | Distribution functions, loop formation probabilities, and force-extension relations in a model for short double-stranded DNA molecules[END_REF], 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 [START_REF] Manghi | Coupling between denaturation and chain conformations in DNA: stretching, bending, torsion and finite size effects[END_REF] and the behavior of stretched DNA molecules [START_REF] Efremov | Transfer-matrix calculations of DNA polymer micromechanics under tension and torque constraints[END_REF] under torsional constraints. Models based on torsional energies alone have also been studied in the context of the statistics of DNA denaturation [START_REF] Fye | Exact method for numerically analyzing a model of local denaturation in superhelically stressed DNA[END_REF][START_REF] Jost | Bubble statistics and positioning in superhelically stressed DNA[END_REF] as well as to model properties of specific DNA sequences [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF].

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 kb) for a wide range of stretching forces and supercoiling levels typical of bacterial genomes in vivo. We next tackle the inference prob-lem 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 [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF] 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.

Theory

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 ( p = K/k B T ), ii) a torsional modulus (C) measured in units of twist persistence length, iii) a winding angle at rest (ψ) of the implicitly embedded helix, iv) the distance (a) between any two consecutive base pairs, thus defining the number of base pairs per p , and v) an effective radius (r e ) reflecting the hard-core simplification of the electrostatic repulsion of the DNA backbone [START_REF] Rybenkov | The effect of ionic conditions on DNA helical repeat, effective diameter and free energy of supercoiling[END_REF]. For B-DNA in physiological conditions ([NaCl] ∼ 100 mM), p typically lies in [40 -60 nm] [62, 63] (in the following, we use 50 nm), C is on the order of 100 nm [START_REF] Strick | Micro-mechanical measurement of the torsional modulus of DNA[END_REF][START_REF] Bryant | Structural transitions and elasticity from torque measurements on DNA[END_REF], while a = 0.34 nm, ψ = 0.6 rad/bp and r e ≈ 2 nm [START_REF] Rybenkov | The effect of ionic conditions on DNA helical repeat, effective diameter and free energy of supercoiling[END_REF].

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 ≤ 10 such that one cylinder corresponds at most to one B-DNA helix. The energy of any conformation C of a self-avoiding linear chain in the presence of a stretching force, f , then reads:

E sRLC (C) = E B (C) + E T (C) -f z, ( 1
)
where z is the extension of the chain along the axis of the force, while

E B (C) = k B T 2 N i=1 p na θ 2 i and E T (C) = k B T 2 N i=1
nC a (φ iψ) 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 θ i is the bending angle associated with the n base pairs of the site and φ i is the winding angle per base pair (Fig. 1A); this means that nφ i is the torsional counterpart of θ 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 A B Figure 1: A. Discrete version of the sRLC: two complete cylinders surrounding a site i are represented along with their local frames ( t, u, v) used to compute the bending angle θ i and the schematically depicted twist angle per bp φ i . B. Extension of the model by allocating a variable s specifying the DNA form of every site. In this case, for symmetry reasons, the length of a cylinder located between sites k -1 and k is equal to n(as k-1 + as k )/2, meaning that cylinders may have different lengths because of the different values of as characterizing the different forms. Here, we depict the case where sites i -1 and i are under the B form (chain in blue), while sites i + 1 and i + 2 are denatured (chain in red). We also indicate by a gray disk the domain wall penalty, J, associated with the structural cost for transiting from one form to the next. the force axis. Altogether, a linear chain made of N bp base pairs thus contains N + 1 sites and N segments (N = N 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 θ i , φ i can be explicitly written as a function of the local frames associated with the surrounding cylinders [START_REF] Carrivain | In Silico Single-Molecule Manipulation of DNA with Rigid Body Dynamics[END_REF] (Fig. 1A). Using Tw ≡ (2π) -1 N i=1 nφ i , this ensures that the linking number, Lk = Tw + 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, 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 [START_REF] Vologodskii | Extension of torsionally stressed DNA by external force[END_REF][START_REF] Lepage | Modeling Bacterial DNA: Simulation of Self-Avoiding Supercoiled Worm-Like Chains Including Structural Transitions of the Helix[END_REF]. These respectively play the role of the fixed surface and of the magnetic bead in single-molecule experiments.

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 [START_REF] Fye | Exact method for numerically analyzing a model of local denaturation in superhelically stressed DNA[END_REF][START_REF] Yan | Localized single-stranded bubble mechanism for cyclization of short double helix DNA[END_REF][START_REF] Wiggins | Exact theory of kinkable elastic polymers[END_REF][START_REF] Ranjith | Distribution functions, loop formation probabilities, and force-extension relations in a model for short double-stranded DNA molecules[END_REF][START_REF] Manghi | Coupling between denaturation and chain conformations in DNA: stretching, bending, torsion and finite size effects[END_REF][START_REF] Efremov | Transfer-matrix calculations of DNA polymer micromechanics under tension and torque constraints[END_REF] and associate a state variable s with every site (s = 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 s ), C s , a s and ψ 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 to be equal to n(a s +a s )/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 [START_REF] Santalucia | A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics[END_REF][START_REF] Huguet | Ritort, Single-molecule derivation of salt dependent base-pair free energies in DNA[END_REF], reflecting the internal deformation of base pairing and base pair stacking. Here, we denote this cost by γ 0 s and set γ 0 B = 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 [START_REF] Peck | Energetics of B-to-Z transition in DNA[END_REF][START_REF] Manghi | Coupling between denaturation and chain conformations in DNA: stretching, bending, torsion and finite size effects[END_REF][START_REF] Oberstrass | Torque Spectroscopy of DNA: Base-Pair Stability, Boundary Effects, Backbending, and Breathing Dynamics[END_REF] (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 C having an arbitrary sequence {s i } i=1..N of DNA forms reads:

E 2sRLC (C) = E γ (C) + E B (C) + E T (C) + E J (C) -f z, (2)
where E γ (C) = n N i=1 γ 0 si is the total free energy formation costs coming from alternative forms,

E B (C) = k B T 2 N i=1 s i nas i θ 2 i is the bending energy, E T (C) = k B T 2 N i=1 nCs i as i (φ i -ψ si ) 2 is the torsional energy and E J (C) = J N -1 i=1 (1 - δ si,si+1
) is the total wall penalty, with δ si,si+1 = 1 if s i = s i+1 , 0 otherwise.

Results and discussion

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 [START_REF] Gebe | Monte Carlo simulations of supercoiling free energies for unknotted and trefoil knotted DNAs[END_REF] and on the statistics of DNA denaturation [START_REF] Fye | Exact method for numerically analyzing a model of local denaturation in superhelically stressed DNA[END_REF][START_REF] Jost | Bubble statistics and positioning in superhelically stressed DNA[END_REF], it is possible to integrate out the torsional degrees of freedom (the φ 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 [START_REF] Vologodskii | Extension of torsionally stressed DNA by external force[END_REF][START_REF] Lepage | Thermodynamics of Long Supercoiled Molecules: Insights from Highly Efficient Monte Carlo Simulations[END_REF] 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):

E ef f γ (C) = n N i=1
γ si where γ si = γ 0 si -

1 2n ln a si C si (3)
E ef f T (C) = k B T 2
N bp χ (C) 2π Tw(C) N bp -ψ(C) 2 + ln( χ (C)) (4)
where, again, N bp = nN is the length of the molecule in base pairs, ψ(C)

≡ N i=1 ψ si N and χ (C) ≡ N i=1 as i Cs i
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

ψ(C) = ψ B and χ (C) = a B /C B for all conformations, such that E ef f T (C) = k B T C B N bp 2a B 2π Tw(C) N bp -ψ B 2
as demonstrated in [START_REF] Gebe | Monte Carlo simulations of supercoiling free energies for unknotted and trefoil knotted DNAs[END_REF].

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 [START_REF] Vologodskii | Simulation of equilibrium and dynamic properties of large DNA molecules, Computational Studies of DNA and RNA[END_REF] in order to include the possibility for mechanical parameters to change at any site of the chain, at any iteration of the simulation (see [START_REF] Lepage | Modeling Bacterial DNA: Simulation of Self-Avoiding Supercoiled Worm-Like Chains Including Structural Transitions of the Helix[END_REF] 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 [START_REF] Lepage | Modeling Bacterial DNA: Simulation of Self-Avoiding Supercoiled Worm-Like Chains Including Structural Transitions of the Helix[END_REF], 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 σ = (Lk -Lk 0 )/Lk 0 . As an example, we report in Fig. 2 both our simulation results and the experimental results of Cees Dekker's lab [START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF] for a 21 kb molecule, for σ ∈ [-0.05, 0.05], f = 0.5 pN and f = 4.5 pN, 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 [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF], 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 nm, D = 15 nm, C D = 10 nm and γ D = 2k B T (see below for discussion of these values) -we intentionally used a D much larger than the 4 nm diameter of the chain because smaller values would lead, in any case, to larger effective persistence lengths [START_REF] Barrat | Persistence length of polyelectrolyte chains[END_REF]. We used J = 10k B T , in accord with previous analyses [START_REF] Peck | Energetics of B-to-Z transition in DNA[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Oberstrass | Torque Spectroscopy of DNA: Base-Pair Stability, Boundary Effects, Backbending, and Breathing Dynamics[END_REF], and considered a discretization level of n = 10, leading to 3.3 cylinders per D and 15 cylinders per B (= 50 nm), a value that both offers reasonable computational times and ensures properties of B-DNA superstructuring to be insensitive to discretization procedures [START_REF] Lepage | Thermodynamics of Long Supercoiled Molecules: Insights from Highly Efficient Monte Carlo Simulations[END_REF][START_REF] Lepage | Modeling Bacterial DNA: Simulation of Self-Avoiding Supercoiled Worm-Like Chains Including Structural Transitions of the Helix[END_REF]. Finally, the quantitative reproduction of the buckling transition at σ > 0 (right part of the curves in Fig. 2) required us to use C = 80 nm at f = 4.5 pN and C = 100 nm at f = 0.5 pN (Fig. S1). This is in accord with the observed systematic discrepancies occurring at low forces between sRLC and experimental extensions [START_REF] Nomidis | Twist-Bend Coupling and the Torsional Response of Double-Stranded DNA[END_REF], a consequence of a coupling between torsional and bending deformations stemming from the asymmetry of the DNA helix [START_REF] Marko | Bending and twisting elasticity of dna[END_REF].

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 0.5 pN) that is symmetric under the transformation σ → -σ and where torsional stress is partially relaxed through the formation of plectonemes, some of them being branched, and ii) a σ-asymmetric high force regime (f 2 pN) where negative torsional stress is partially relaxed by forming denaturation bubbles. For intermediate forces at σ < 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 [START_REF] Salerno | Single-Molecule Study of the DNA Denaturation Phase Transition in the Force-Torsion Space[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF], 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 [START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. 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 [START_REF] Laio | Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science[END_REF], 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 [START_REF] Matek | Plectoneme tip bubbles: Coupled denaturation and writhing in supercoiled DNA[END_REF]. Here, we observe that these tip-bubble conformations, as coined in [START_REF] Matek | Plectoneme tip bubbles: Coupled denaturation and writhing in supercoiled DNA[END_REF], 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 [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF][START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], while previous works have led to values of the torsional modulus that may differ by more than 20-fold [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. 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.

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 ( D , C D , a D , ψ D and γ 0 D ) 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 pN and 2 pN cannot thus be currently exploited to further constrain parameters as previously done in the context of phenomenological models [START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF].

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 Lk = 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,σ ({λ i }), as a function of the fraction of each of the possible alternative forms at a given stretching force f and supercoiling density σ. In the following, for simplicity we discuss the case of two forms, with in mind B-DNA plus a certain fraction λ of an alternative form, hereafter denoted as X-DNA -the generalization to more than two forms is straightforward (see Supplementary methods).

Calling J (λ) the free energy associated with the multiple possibilities to distribute the λ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 χ s = a s /C s and χ /ψ(λ) = (1λ) χ B /ψ B + λ χ s /ψ s , we obtain (Supplementary methods):

F f,σ (λ) = λN bp γ X + J (λ)
free energy formation of alternative domains

+ g X,f (λN bp a X ) + g B,f ((1 -λ)N bp a B ) bending+stretching + k B T 2
N bp χ (λ) 1 + σ - ψ(λ) ψ B 2 + ln ( χ (λ)) torsion (5)
Using this free energy profile, equilibrium properties of the 2sRLC model in the zero writhe approximation can be computed by performing numerical integration over λ 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 * = -∂ f F f,σ (λ) λ , the torque Γ * = (2πTw 0 ) -1 ∂ σ F f,σ (λ) λ (see Supplementary methods for derivation and further details), and the fraction of alternative form, λ * = λ λ , where

• λ ≡ dλ • exp [-F f,σ (λ)/k B T ] dλ exp [-F f,σ (λ)/k B T ].
We also compute the equilibrium number of domains, X * , using a free energy surface in the (λ, X) space (Eq. 6 in Supplementary methods).

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 X = 15 nm, a X = 0.54 nm, ψ X = 0, C X = 10 nm and γ X = 1.6kT , having in mind a form close to D-DNA; just as in the analysis of Fig. 2, the relatively high value of X compared to previous estimations for D-DNA (3 -4 nm [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF][START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF]) 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 [START_REF] Lepage | Modeling Bacterial DNA: Simulation of Self-Avoiding Supercoiled Worm-Like Chains Including Structural Transitions of the Helix[END_REF] and to avoid considering (longer) effective persistence lengths in our analytical 3: Comparison between numerical simulations, the zero writhe approximation and Marko's phenomenological approach for the equilibrium value of different quantities of interest as a function of the negative supercoiling: the fraction λ * of denatured alternative form, the number X * of alternative domains, the torque Γ * and the extension of the molecule (relative to its contour length under the single B form). Error bars represent one standard deviation due to thermal fluctuations; in particular, this reveals that torque fluctuations are large. For readability, the green points (2 pN) and red points (10 pN) have respectively been slightly shifted left and right with respect to the exact values of σ: 0, -0.1, -0.2, -0.3, -0.4, -0.5 (-0.015, -0.02, -0.025, -0.03 in the inset of X * ). The transition between fully stretched conformations and plectonemic conformations (where the approximation does not hold anymore) occurs between 1 and 2 pN. For all panels: the apparent absence of error bars implies that they are smaller than symbols.

approximations [START_REF] Barrat | Persistence length of polyelectrolyte chains[END_REF]. We also used a fine discretization level of n = 4 to get a large number of cylinders per D ( 7) so that to prevent any discrepancy coming from the discretization procedure of the rod-like chain.

Regarding λ * and the extension, we observe an excellent agreement between simulation results and the zero-writhe approximation for any force above 2 pN (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 σ -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 σ = -0.5, and for forces down to the onset of plectonemic superstructures. The strong deviation observed at 1 pN 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 [START_REF] Seol | Elasticity of Short DNA Molecules: Theory and Experiment for Contour Lengths of 0.6-7µm[END_REF] 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 [START_REF] Bouchiat | Estimating the persistence length of a worm-like chain molecule from force-extension measurements[END_REF], 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 [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF]. 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 λ of X-DNA and the total supercoiling density (σ) 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 σ. 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 [START_REF] Marko | Torque and dynamics of linking number relaxation in stretched supercoiled DNA[END_REF].

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 [START_REF] Santalucia | A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics[END_REF] and singlemolecule [START_REF] Huguet | Ritort, Single-molecule derivation of salt dependent base-pair free energies in DNA[END_REF] 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 φ D = 0 [START_REF] Kahn | Detection of localized DNA flexibility[END_REF]. A zero-writhe phenomenological modeling of these experiments have then led to a D = 0.54 nm, D 2 to 3 nm and C D = 1 nm [START_REF] Kahn | Detection of localized DNA flexibility[END_REF]. In contrast, investigation of the co-existence of D-DNA with B-DNA at forces 1 pN and |σ| 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 nm [START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. Regarding left-handed DNA forms, static and dynamic laser light scattering experiments of Z-DNA solutions led to Z 200 nm [START_REF] Thomas | Chain flexibility and hydrodynamics of the B and Z forms of poly(dG-dC).poly(dG-dC)[END_REF], whereas a best fitting procedure using Marko's phenomenological model revealed a flexible L-DNA form, with a L = 0.48 nm and L = 3 nm [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], and a torsional modulus C L lying between 10 and 20 nm. Similar values of C L were found by investigating the collapse of pure L-DNA molecules [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF]. Finally, location of the buckling transition of such "pure" L-DNA molecules at high forces led to ψ L = -0.5 rad/bp [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF].

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 [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], to estimate parameters of D-DNA and L-DNA using σ-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.

Simulated data: assessing parameter inference in an underdetermination

context To test our capacity to recover original parameters (hereafter indicated by a hat, e.g., âX ) from simulated σ-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 [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], 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 X (or equivalently K X ) have a more systematic impact on σ-extension curves. More precisely, using a large random set of parameters, the various σ-extension curves appear to be sorted according to the values of both a X and 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 ψ X that are distant from ψ 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 , 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 (â X , ˆ 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 X , the exact value being irrelevant); 2) solutions corresponding to more disordered domains (small X ) with longer double helices (larger a X ). Most importantly, an analysis of solutions in the vicinity of (â X , ˆ X ) corroborates the existence of a large spectrum of possible values for J, γ X , ψ 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 âX and ˆ X , too.

In principle, the value of ĈX could be addressed using torque measurements [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF], yet only if one has access to λ(σ) (see Supplementary methods). For âX and ˆ 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 γ X leads to the same underdeter-Figure 4: Test of the possibility to recover parameters of an alternative X-DNA co-existing with B-DNA, using simulated extension curves (red curves in panel A). To this end, we first generated σ-extension curves using the zero writhe approximation model of Eq. 5 with a specifically predefined X-DNA (see main text for parameters). We next performed an inference procedure using the same generative model by producing 10 000 different curves corresponding to a random combination of the 6 parameters of X-DNA. A: σ-extension curves at f = 2 pN where the color of each curve is given by the value of one of the parameters: a X (left), X (center) and γ X (right). See Fig. S6 for the other parameters. We can see that, compared to γ X , a X and X have a well-defined, strong influence on extension curves. In these cases, the curves indeed appear to be sorted as a function of the values of the parameters. B: Each point corresponds to a random set of parameter values, with x-and y-axes corresponding to the values of a X and X , respectively. The opacity is inversely correlated (using an arbitrarily scale) to the root-mean-square deviation (RMSD) of the curves in A with respect to the red one, so that only the best fits are visible. We observe a crest of small RMSD values that go through the original parameters (â X , ˆ X ) (black cross). C: At a higher force (12 pN), the front is much less visible and the original set of solutions becomes isolated among a broad cloud of RMSD values (in blue). The crest becomes clearer again by setting ψ X = ψX (= 0) and C X = ĈX (= 10 nm) (in red). The intersection between crests obtained at different forces then allows for an accurate estimation of âX and ˆ X (Fig. S7). mination problem. In contrast, fixing either ψ X or C X makes a crest of locally suboptimal solutions re-appear. However, setting ψ X to the right value ψX = 0 results in a crest that is off the right solution (â X , ˆ X ) (Fig. S9), a shift that disappears only by setting C X = ĈX (red points in Fig. 4C).

Overall, these results show that inference of âX and ˆ X from extension curves requires to have some prior knowledge about the nature of the alternative form, more particularly about its torsional parameters ( ψX and ĈX ). This may be explained by the fact that the (effective) torsional energy is quadratic in σ, such that any parameter that affects this energy is expected to have a large impact on any curve plotted as a function of σ. When the values of ψX and ĈX are known, a minimization over all forces then lead to a fair estimation of both âX and ˆ 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).

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

To estimate â and ˆ p for both D-DNA and L-DNA, we use the data of Sheinin et al. [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF] obtained for 2.2 kb long molecules, which consist of extension curves measured over a wide range of values of σ (from -2.3 to +0.1) and f (from 2.5 to 36 pN). In accord with the findings of [START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF], we further consider that L-DNA occurs for σ -1 and, hence, estimate D-DNA and L-DNA parameters using an optimization procedure of the extension curves for σ 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, ψ D = 0 and ψ L = -0.5 rad/bp and we test, on the other hand, C D = 1 nm as discussed in [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], C D = 28 nm as proposed in [START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF] and C L = 10 nm in accord with previous estimations [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF].

Similar to the case of simulated data, optimizations over forces from 2.5 to 12 pN 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 nm leads to a D > 0.6 nm and D close to 1 nm, which are different, although being qualitatively similar, from previous estimations according to which a D = 0.54 nm and D ∈ [2-3 nm] (see [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF] and references therein). When C D = 28 nm, we find a D ∈ [0.45 -0.5 nm] and D ∈ [2 -5 nm]. 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 J (λ) 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 âD and ˆ 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 ∈ [0.5 -0.6 nm], which is slightly larger than previous estimations [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF], and L ∈ [1 -2 nm], which is slightly smaller. In all cases, the small value of L corroborates that L-DNA is not a pure left-handed structural form but made of denatured DNA [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF]. [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF]. Each point represents the cumulative (over several forces ) RMSD resulting from a best fit procedure of experimental extension curves using numerical curves of the 2sRLC model generated under the zero-writhe approximation (see Eq. 5) and for a given set of parameters. Points are all the darker that the cumulative RMSD is small (arbitrary scale). For each set of parameters, the cumulative RMSD was computed for σ varying in [-0.6, -0.1] and over 6 different forces ranging from 2.5 to 12 pN. To mitigate underdetermination problems, we set ψ D = 0 and tested two values of C D , which are representative of values discussed in the literature [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. For C D = 1 nm (red points), we find an island of optimal solutions characterized by a D > 0. To summarize, our results question the compatibility of previously used a D = 0.54 nm 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.

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 [START_REF] Cong | Recent progress on the mechanics of sharply bent DNA[END_REF] and the superstructuring properties of supercoiled B-DNA molecules [START_REF] Strick | Stretching of macromolecules and proteins[END_REF][START_REF] Vologodskii | Biophysics of DNA[END_REF]. 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 ATtracks and CpG-tracks to respectively favor DNA denaturation and left-handed forms [START_REF] Mirkin | Discovery of alternative DNA structures: a heroic decade (1979-1989)[END_REF], can be easily realized by using sequence dependent free energy formation costs of alternative forms (the γ si 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 γ 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 [START_REF] Matek | Plectoneme tip bubbles: Coupled denaturation and writhing in supercoiled DNA[END_REF], we find that for a force f = 1 pN and a supercoiling density σ = -0.06, the site i AT is systematically denatured and serves as both nucleation and anchor points for a plectoneme (Fig. S11).

Conclusion

More elaborate 2sRLC models can be contemplated by e.g. including sequence effects at the level of DNA bending [START_REF] Geggier | Sequence dependence of DNA bending rigidity[END_REF][START_REF] Mitchell | Maddocks, Sequence-Dependent Persistence Lengths of DNA[END_REF]. In particular, by properly distinguishing intrinsic curvature from stiffness [START_REF] Mitchell | Maddocks, Sequence-Dependent Persistence Lengths of DNA[END_REF], one can expect to quantitatively address the structuring properties of bacterial sequences, including the possible pinning of pure B-DNA plectonemes [START_REF] Kim | DNA sequence encodes the position of DNA supercoils[END_REF]. 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 [START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. 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 [START_REF] Efremov | Transfer-matrix calculations of DNA polymer micromechanics under tension and torque constraints[END_REF]. 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 [START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF].

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 [START_REF] Oberstrass | Torque measurements reveal sequence-specific cooperative transitions in supercoiled DNA[END_REF][START_REF] Oberstrass | Torque Spectroscopy of DNA: Base-Pair Stability, Boundary Effects, Backbending, and Breathing Dynamics[END_REF], 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 (λ 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. [START_REF] Oberstrass | Torque Spectroscopy of DNA: Base-Pair Stability, Boundary Effects, Backbending, and Breathing Dynamics[END_REF]) and, hence, should constrain further model predictions -just as for the co-existence of plectonemic conformations and strecthed conformations with alternative forms (Figs. 2 andS2), 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 [START_REF] Vlijm | Experimental phase diagram of negatively supercoiled DNA measured by magnetic tweezers and fluorescence[END_REF], might offer solutions to compute the fraction λ.

Using Tw(C) = (2π) -1 i nφ i and defining t(C) = Tw Lk (C)/N bp = Tw Lk (C)/nN , Eq. 1 can be rewritten as:

P (C) = Z -1 exp -βE T (C) - i nC s i 2a s i (φ i -ψ s i ) 2 δ i φ i -2πtN (2)
The integration over the φ i 's then results in evaluating the following term (for clarity, we drop the C's):

I = k dφ k exp - i nC s i 2a s i (φ i -ψ s i ) 2 δ i φ i -2πtN = (using the Fourier transform of δ) dx k dφ k exp - nC s k 2a s k (φ k -ψ s k ) 2 -2πix(2πt -φ k ) = dx k dφ k exp - nC s k 2a s k φ k -ψ s k - 2πia s k nC s k x 2 exp - 2π 2 a s k nC s k x 2 -2πix(2πt -ψ s k ) = k dφ k exp - nC s k 2a s k φ k -ψ s k - 2πia s k nC s k x 2 dx exp - j 2π 2 a s k j nC s j x 2 + 2πix(2πt -ψ s j ) (3)
Defining the conformation-dependent average torsional susceptibility per base pair χ (C) ≡ N i=1 as i Cs i

N , the conformation-dependent average twist angle at rest ψ(C) ≡ N i=1 ψ s i N , we obtain (still dropping the C's):

I = (2π) N/2 k a s k nC s k dx exp - 2π 2 N χ n x 2 + i n π χ (2πt -ψ)x = (2π) N/2 k a s k nC s k dx exp - 2π 2 N χ n x + i n 2π χ (2πt -ψ) 2 + n 2 4π 2 χ 2 (2πt -ψ) 2 = (2π) N/2 k a s k nC s k exp - nN 2 χ (2πt -ψ) 2 dx exp - 2π 2 N χ n x + i n 2π χ (2πt -ψ) 2 = (2π) (N -1)/2 N 1/2 k a s k χ C s k exp - nN 2 χ (2πt -ψ) 2 = (2π) (N -1)/2 N 1/2 exp - 1 2 N bp 2 χ (2πt -ψ) 2 + ln( χ ) - i ln a s i C s i , (4)
which, after redistribution of the terms ln( χ ) and ln Using Eq. 6, one can then compute the equilibrium number of domains in the zero writhe approximation by numerically determining:

X * = dXdλ X exp [-S f,σ (λ, X)/k B T ] dXdλ exp [-S f,σ (λ, X)/k B T ] (10)
In addition, using Eq. 7 and

• λ ≡ dλ • exp [-F f,σ (λ)/k B T ] dλ exp [-F f,σ (λ)/k B T ],
one can derive the other equilibrium quantities as follows (see below for torques):

λ * = λ λ , (11)
z * = -∂ f F f,σ (λ) λ = λ * X + (1 -λ * ) B (12)
with s the WLC extension associated with the form s of DNA, which depends only on a s and s [1].

Note finally that numerical integration is performed by circumventing numerical problems coming from the generation of large combinatorial quantities when λN becomes large (as in the case e.g. of Eq. 9). Recalling that λ = K/N where K is the number of sites occupied by the alternative form, the integration thus consists of a discrete sum over restricted values of K and X for which Boltzmann weights are maximal. To this end, we first localise the saddle-point of S f,σ (λ, X) and then perform the sum over values of K and X lying in the rectangle that contains the saddle-point and outside of which relative Boltzmann weights with respect to the saddle-point are smaller than 10 -4 . We further checked that our result did not depend on the choice of this cut-off.

C. Generalization to an arbitrarily number of alternative forms

In the situation of k possible DNA forms (one B-DNA plus k -1 alternative forms), denoting the vector of fraction of each of these forms by λ ≡ {λ i } i=1..k , with k i=1 λ i = 1, the free energy landscape becomes:

F f,σ (λ) = N bp k i=1 λ i γ i + I(λ) + k i=1 g i,f (λ i N bp a i ) + k B T 2   ψ 2 B N bp k i=1 λ i χ i 1 + σ - k i=1 λ i ψ i ψ B 2 + ln k i=1 λ i χ i   (13)
where here the index i always indicates the different forms and where I(λ) indicates the free energy associated with the multiple possibilities to distribute the different alternative sites into distinct domains.

D. Torques

Measurement in numerical simulations

To measure the torque exerted by a molecule in our simulations, we use the same method as in our previous work [2] where we adapted a commonly used experimental method [3,4] to the case of rod-like chain models with global torsional energy. It consists in attaching the ends of the molecule to a reservoir of twists with which this can exchange helix turns. The reservoir, also called a magnetic trap, is characterized by a stiffness k R such that any variation of the twist within it, δTw R , with respect to a reference twist (see below) is associated with an energy cost 2π 2 k B T k R δTw 2 R . Similarly to the measurements in magnetic torque tweezers, the equilibrium torque exerted by the molecule is then given by Γ * = -2πk R (Tw * -Tw * 0 ), where Tw * and Tw * 0 are respectively the equilibrium twists of the supercoiled molecule and of the corresponding relaxed molecule, both measured in the presence of the magnetic trap.

Note here, that the linking number in the reservoir is further calibrated such that Tw * 0 is equal to the number of helices of the molecule at rest, Tw 0 .

Expression in the zero writhe approximation and access to the torsional modulus

The equilibrium torque exerted by the molecule is defined by Γ * = ∂ Θ F f,σ (λ) λ where Θ is the rotation angle of the tip of the molecule around the stretching direction (the direction perpendicular to the magnetic field used to add helical turns to the molecule). In magnetic tweezers experiments, Θ is a control parameter related to the linking number: Θ = 2πLk.

As a consequence, Γ * = (2π) -1 ∂ Lk F f,σ (λ) λ , which becomes Γ * = (2πTw 0 ) -1 ∂ σ F f,σ (λ) λ in the approximation of zero writhe, where Lk = Tw and σ = (Tw -Tw 0 )/Tw 0 . Using Eq. 7, we thus obtain for the case of the co-existence of X-and B-DNA:

Γ * = k B T ψ B 1 + σ -ψ(λ) ψ B χ (λ) λ (14)

Figure 2 :

 2 Figure 2: Extension of a 21 kb molecule as a function of σ, for different forces, obtained experimentally (in red and blue, data from [27]) and numerically using a 2sRLC model that includes the possibility to have B-DNA plus an alternative DNA form along the chain (see main text for values of the microscopic parameters). Error bars and colored areas represent one standard deviation due to thermal fluctuations. At low forces, the curves are symmetrical with respect to σ, with the formation of branched plectonemic structures at large values of |σ| as indicated by the typical conformations observed in the simulations (leftmost and rightmost panels). At high forces, the buckling transition toward the plectonemic regime occurs only for positive supercoiling. For negative supercoiling, denaturation bubbles appear instead (upper left panel, in red). Upper panel: corresponding estimation of the equilibrium fraction λ * of the alternative form as a function of σ.

Figure

  Figure3: Comparison between numerical simulations, the zero writhe approximation and Marko's phenomenological approach for the equilibrium value of different quantities of interest as a function of the negative supercoiling: the fraction λ * of denatured alternative form, the number X * of alternative domains, the torque Γ * and the extension of the molecule (relative to its contour length under the single B form). Error bars represent one standard deviation due to thermal fluctuations; in particular, this reveals that torque fluctuations are large. For readability, the green points (2 pN) and red points (10 pN) have respectively been slightly shifted left and right with respect to the exact values of σ: 0, -0.1, -0.2, -0.3, -0.4, -0.5 (-0.015, -0.02, -0.025, -0.03 in the inset of X * ). The transition between fully stretched conformations and plectonemic conformations (where the approximation does not hold anymore) occurs between 1 and 2 pN. For all panels: the apparent absence of error bars implies that they are smaller than symbols.

Figure 5 :

 5 Figure5: A: Inference of D-DNA parameters using experimental data from[START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF]. Each point represents the cumulative (over several forces ) RMSD resulting from a best fit procedure of experimental extension curves using numerical curves of the 2sRLC model generated under the zero-writhe approximation (see Eq. 5) and for a given set of parameters. Points are all the darker that the cumulative RMSD is small (arbitrary scale). For each set of parameters, the cumulative RMSD was computed for σ varying in [-0.6, -0.1] and over 6 different forces ranging from 2.5 to 12 pN. To mitigate underdetermination problems, we set ψ D = 0 and tested two values of C D , which are representative of values discussed in the literature[START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. For C D = 1 nm (red points), we find an island of optimal solutions characterized by a D > 0.6 nm and D close to 1 nm. For C D = 28 nm, we find a D ∈ [0.45-0.5 nm] and D ∈ [2-5 nm]. The dashed horizontal line in the figure indicates D = 2 nm. B: Same analysis but for L-DNA. Here, we set ψ L = -0.5 rad/bp and test C L = 10 nm [25, 26]. We find a L ∈ [0.5 -0.6 nm] and L ∈ [1 -2 nm].

  Figure5: A: Inference of D-DNA parameters using experimental data from[START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF]. Each point represents the cumulative (over several forces ) RMSD resulting from a best fit procedure of experimental extension curves using numerical curves of the 2sRLC model generated under the zero-writhe approximation (see Eq. 5) and for a given set of parameters. Points are all the darker that the cumulative RMSD is small (arbitrary scale). For each set of parameters, the cumulative RMSD was computed for σ varying in [-0.6, -0.1] and over 6 different forces ranging from 2.5 to 12 pN. To mitigate underdetermination problems, we set ψ D = 0 and tested two values of C D , which are representative of values discussed in the literature[START_REF] Sheinin | Underwound DNA under Tension: Structure, Elasticity, and Sequence-Dependent Behaviors[END_REF][START_REF] Meng | Coexistence of Twisted, Plectonemic, and Melted DNA in Small Topological Domains[END_REF]. For C D = 1 nm (red points), we find an island of optimal solutions characterized by a D > 0.6 nm and D close to 1 nm. For C D = 28 nm, we find a D ∈ [0.45-0.5 nm] and D ∈ [2-5 nm]. The dashed horizontal line in the figure indicates D = 2 nm. B: Same analysis but for L-DNA. Here, we set ψ L = -0.5 rad/bp and test C L = 10 nm [25, 26]. We find a L ∈ [0.5 -0.6 nm] and L ∈ [1 -2 nm].

  as i Cs i in the torsional energy and the free energy of formation, respectively, leads to the effective energies of the main text (Eqs. 3 and 4). Note also that in Eq. 4 of the main text, we have used the fact that Tw(C) = Tw Lk (C) in the simulated polymer model. the number of ways to distribute Y sites in X different domains among N possible sites.

  FIG. S2. Hysteresis cycles around the plectonemic-denatured transition. Simulations of a 2.4 kb molecule were performed (σ = -0.05) during which the force was continuously increased (as indicated by the lower left arrow) and decreased (upper right arrow arrow) at constant speed from 0.3 to 1.8 pN. For each of the three speeds reported here (from the fastest to the slowest: red, green and blue curves), 12 cycles were performed, revealing in each case a strong hysteresis pattern in the force-extension diagram.

  FIG.S9. Setting only ψ X = ψX (= 0), at high force (12 pN) we find a crest of sub-optimal fits that go off the original parameters (â X , ˆ X ) (black cross).

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.

Supplemental material: A polymer model of bacterial supercoiled DNA including structural transitions of the double helix In the following, parameters are uniform along the DNA molecule. Sequence effects will be treated in subsequent works.

A. Derivation of the effective model with a global torsional energy

We start by writing down the equilibrium probability of a configuration C with DNA-form profile {s i } i=1..N in the context of an imposed linking number Lk. To this end, we divide the conformational energy into two terms: the torsional energy

that depends on the torsional angles {φ i } i=1..N and the rest, E T (C), which does not. In this context, denoting β = 1/k B T , we have:

with Z the partition function such that C P (C) = 1. The δ(x) function, which is equal to 1 if x = 0 and 0 otherwise, imposes the conservation of Lk, that is, it only retains conformations for which the total twist of the molecule, Tw(C), is equal to Tw Lk (C) ≡ Lk -Wr(C).

B. Free energy profiles in the zero writhe approximation

In the case of the co-existence of two forms, here B-DNA and X-DNA (any alternative DNA form), using Eqs. 2, 3 and 4 of the main text, the effective energy, U (λ, X), of a conformation where X-DNA occupies λN sites (corresponding to λN bp base pairs) distributed in X domains reads:

where E s B indicates the total bending energy of the form s and where we have used the notations χ s = a s /C s and χ /ψ(λ) = (1λ) χ B /ψ B + λ χ s /ψ s -note that ψ(λ) represents the average twist angle at rest of the molecule with a fraction λ of alternative forms. Next, considering that the supercoiling density, σ, is totally converted into twist in the approximation of zero writhe, we can replace 2π Tw N bp by (1 + σ)ψ B . Finally, integrating the corresponding partition function over the bending angles of each site leads to a free energy where the terms E s Bf z s are replaced by g s,f (X), the free energy associated with the WLC at force f with bending modulus K s and contour length X. The resulting bi-dimensional free energy surface, S f,σ (λ, X) therefore reads:

The unidimensional free energy profile F f,σ (λ) (Eq. 5 in the main text) can be obtained by summing (with Boltzmann weights) over the values of X, leading to:

with:

where we have defined where the cylinder at the center (number 50) has no denaturation penalty (γ D,i=50 = 0). The denaturation bubble (red) appears and stays at the center during the whole simulation.