# Observation of surface superconductivity in a three-dimensional Dirac
material

Qi Liu  1,† Peng-Jie Guo  2,† Xiao-Yu Yue  1 Zhe-Kai Yi  1 Qing-Xin Dong  5
Hui Liang  1 Dan-Dan Wu  1 Yan Sun  1 Qiu-Ju Li  4 Wen-Liang Zhu  6 Tian-Long
Xia  2,‡ Xue-Feng Sun3,1,7,§ and Yi-Yan Wang1,∗

###### Abstract

Superconductivity becomes more interesting when it encounters dimensional
constraint or topology, because it is of importance for exploring exotic
quantum phenomena or developing superconducting electronics. Here we report
the coexistence of naturally formed surface superconducting state and three-
dimensional topological Dirac state in single crystals of BaMg2Bi2. The
electronic structure obtained from the first-principles calculations
demonstrates that BaMg2Bi2 is an ideal Dirac material, in which the Dirac
point is very close to the Fermi level and no other energy band crosses the
Fermi level. Superconductivity up to 4.77 K can be observed under ambient
pressure in the measurements of resistivity. The angle dependent
magnetoresistance reveals the two-dimensional characteristic of
superconductivity, indicating that superconductivity occurs on the surface of
the sample and is absent in the bulk state. Our study not only provides
BaMg2Bi2 as a suitable platform to study the interplay between
superconductivity and topological Dirac state, but also indicates that MgBi-
based materials may be a promising system for exploring new superconductors.

Institute of Physical Science and Information Technology, Anhui University,
Hefei, Anhui 230601, China

Department of Physics and Beijing Key Laboratory of Opto-electronic Functional
Materials & Micro-nano Devices, Renmin University of China, Beijing 100872,
China

Department of Physics and Key Laboratory of Strongly-Coupled Quantum Matter
Physics (CAS), University of Science and Technology of China, Hefei, Anhui
230026, China

School of Physics & Materials Science, Anhui University, Hefei, Anhui 230601,
China

Institute of Physics and Beijing National Laboratory for Condensed Matter
Physics, Chinese Academy of Sciences, Beijing 100190, China

School of Physics and Information Technology, Shaanxi Normal University, Xi’an
710062, China

Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing, Jiangsu 210093, China

†These authors contributed equally to this work.

‡Corresponding authors<EMAIL_ADDRESS>
§Corresponding authors<EMAIL_ADDRESS>
∗Corresponding authors<EMAIL_ADDRESS>
The study of two-dimensional (2D) superconductivity provides insight into a
variety of quantum phenomena in condensed matter physics and material
science[1]. The 2D superconductivity can be realized in systems such as
heterogeneous interfaces or atomically ultrathin materials[1]. The 2D
interfacial superconductivity has been observed in LaAlO3/SrTiO3[2],
FeSe/SrTiO3[3, 4] and EuO/KTaO3(111)[5] _etc_. Few-layer stanene and 1 _T_
d-MoTe2 also exhibit superconductivity at low temperature[6, 7]. In addition,
superconductivity up to 14.9 K was discovered in thin film of infinite-layer
nickelate Nd0.8Sr0.2NiO2[8], indicating the possibility of a family of
nickelate superconductors. Recently, surface superconductivity can also be
found in some topological materials[9], such as Dirac semimetal KZnBi[10],
Cd3As2 nanoplate (proximity-induced)[11] and type II Weyl semimetal
TaIrTe4[12], which provide more platforms for studying 2D superconductivity.

Searching for superconductivity in topological materials has attracted a lot
of interest for the possibility of obtaining novel topological phases, for
example, the topological superconductivity (TSC)[13, 14]. Topological
superconductors host superconducting bulk state and gapless Majorana surface
state[15, 16]. The Majorana fermions exist in the superconducting gap and have
the potential to be applied to topological quantum computation. Recently,
Zhang _et al._ reported the observation of TSC on the surface states of
FeTe0.55Se0.45[17]. Liu _et al._ demonstrated the Majorana zero mode arising
from the chiral topological surface state of (Li0.84Fe0.16)OHFeSe[18]. The TSC
can also be induced by the proximity effect between an _s_ -wave
superconductor and surface Dirac fermions of topological insulator[19].
Experimentally, many attempts have been made to induce superconductivity on
the surface of topological materials, including the fabrication of
heterostructure[20, 21, 22, 23, 24] and tip/proximity-induced
superconductivity[11, 25, 26, 27] _etc_.

In this work we grow the single crystals of BaMg2Bi2 and study the transport
properties and electronic structure. The first-principles calculations show
that BaMg2Bi2 is a three-dimensional (3D) Dirac material, in which the Dirac
point exists in the direction of $\Gamma$-A and is protected by the symmetry.
The Dirac point is very close to the Fermi level. In addition, only a small
hole-type Fermi surface carrying the Dirac fermions crosses the Fermi level in
the band structure, indicating that BaMg2Bi2 is an ideal material to explore
extraordinary properties associated with Dirac fermions. Interestingly,
superconductivity emerges below 4.77 K in the measurements of resistivity. The
polar plots of the angle dependent magnetoresistance (MR) show a two-fold
symmetry in the superconducting state, suggesting the 2D characteristic of the
observed superconductivity. Our findings demonstrate the coexistence of
surface superconductivity and 3D topological Dirac state in BaMg2Bi2.

Figure 1a illustrates the crystal structure of BaMg2Bi2, which consists of
alternating Ba atoms and MgBi layers. BaMg2Bi2 crystalizes in the hexagonal
CaAl2Si2-type structure with the space group of $P\bar{3}m1$ (No. 164). The
crystal structure has C3 rotational symmetry, as shown in the view from _c_
-axis direction (Fig. 1b). The first-principles calculations have been
employed to study the electronic structure of BaMg2Bi2. When the spin-orbital
coupling (SOC) effect is ignored, the band structure exhibits a gap around the
$\Gamma$ point (Fig. S7 in the Supplementary Information). As shown in Figs.
1c and 1d, when the SOC effect is included, a Dirac point appears in the
$\Gamma$-A direction, which is protected by the C3 symmetry. It should be
noted that the Dirac point is very close to the Fermi level and locates only
about 8.5 meV above the Fermi level. Except for the Dirac fermion, there is no
other energy band near the Fermi level. From the band structure, Dirac fermion
will dominate the transport properties of BaMg2Bi2. The calculated band
structure is consistent with the previous report on angle-resolved
photoemission spectroscopy (ARPES) experiments[28], revealing the 3D
topological Dirac state in BaMg2Bi2.

Figure 2a presents the temperature dependent resistivity $\rho$(_T_) of
BaMg2Bi2. The applied current is 100 $\mu$A. At first, the resistivity
decreases with the decrease of temperature, exhibiting the metallic behavior.
After the minimum at 160 K, an upturn appears in the resistivity curve with
decreasing temperature, similar to a metal-insulator transition. Surprisingly,
a sharp resistivity drop shows up at low temperature, indicating the
transition to superconducting state. The superconducting transition
temperature determined by linear epitaxy of normal curve and falling curve is
4.77 K (inset in Fig. 2a). The observed superconducting state has been further
characterized by changing magnetic field and current intensity. Figure 2b
shows the temperature dependent resistivity of BaMg2Bi2 under various magnetic
fields. The applied field is perpendicular to the current direction.
Superconductivity is gradually suppressed with increasing the field. At 8 T
and 9 T, there is only a small drop in resistivity. The temperature dependent
resistivity with different currents is shown in Fig. 2c. The superconducting
transition temperature shifts toward lower temperatures with increasing
current. When the current reaches 4 mA, the superconductivity is almost
completely suppressed.

Figures 3a and 3b display the temperature dependence of critical field _H_ c2
and critical current _I_ c extracted from Figs. 2b and 2c, respectively. As
shown by the fit (red solid line) in Fig. 3a, the variation of critical field
can be well described by the empirical formula _H_ c2 = _H_ c2(0)[1 – (_T_
/_T_ c)2]. The obtained upper critical field at 0 K is 8.7 T. With the
increase of current, _T_ c decreases rapidly. The behavior is similar to the
interfacial superconductivity in TaAs/Ag heterostructure and LaAlO3/SrTiO3
system.

Since BaMg2Bi2 exhibits a superconducting transition in the measurement of
resistivity, the Meissner effect and specific heat jump are expected. However,
after testing on several samples, it turns out that although superconducting
transition can be observed in resistivity, the Meissner effect and a jump in
heat capacity are absent (Figs. S2 and S3 in the Supplementary Information).
This implies the possibility of surface superconductivity in the samples of
BaMg2Bi2. Superconductivity only occurs on the surface of the sample, but the
bulk state accounts for a large proportion in the measurements of magnetic
susceptibility and specific heat.

In order to further understand the unusual superconductivity in BaMg2Bi2, the
measurements of angle dependent resistivity under magnetic field are
performed. The schematic diagram of the measurements are plotted in Figs. 4a
and 4d. The current is along the hexagonal edge of crystal and the magnetic
field rotates in the _yz_ -plane and _xz_ -plane, respectively. In the former
configuration, as shown in Figs. 4b and 4c, the polar plots of the angle
dependent MR exhibit two-fold symmetry with a period of $\pi$. The maximum and
minimum values of resistivity appear at 0∘, 180∘ (_B_ $\bot$ _ab_ plane) and
90∘, 270∘ (_B_ // _ab_ plane), respectively. In the latter configuration, the
patterns in Figs. 4e and 4f show a slight distortion, but the two-fold
symmetry remains. It should be noted that the two-fold symmetry only exists in
the superconducting state. Once the superconducting state is gradually
suppressed with the increase of magnetic field or temperature, the two-fold
symmetry will be weakened and disappear. In the nonsuperconducting state, the
angle dependent MR shows isotropic characteristic, which comes from the 3D
Fermi surface of the bulk state. In general, the two-fold symmetry means that
the conductive layer is 2D. In materials, the Lorentz force affects the
carrier momentum component in the plane perpendicular to the magnetic field.
For a 2D conductive layer, the carrier only responds to the magnetic field
component _B_ $\mid$cos$\theta$$\mid$, resulting in the two-fold symmetry of
the polar plots. In addition, the observed two-fold symmetry cannot be derived
from the bulk state, because the bulk state is a small 3D Fermi surface, as
demonstrated by the first-principles calculations and previous ARPES
experiment. Thus, the superconductivity in BaMg2Bi2 is considered to be a 2D
superconductivity that occurs on the surface of the sample.

In fact, the properties of BaMg2Bi2 are similar to KZnBi, which also hosts the
surface superconductivity (0.85 K) and Dirac state[10]. In both samples, the
superconductivity naturally forms on the surface without external stimuli
while the bulk state remains nonsuperconducting. Since Dirac fermions exist in
the bulk state, it is worth exploring the relationship or interplay between
the superconductivity and topological state. Further researches (for example,
ARPES experiment in superconducting state) are needed to clarify the remaining
problems. On the material side, BaMg2Bi2 possesses the MgBi layers that are
important for transport properties and the formation of Dirac state. The MgBi-
based materials may provide more opportunities for exploring new
superconductors.

In summary, we report the observation of surface superconductivity up to 4.77
K under ambient pressure in the single crystals of BaMg2Bi2, which hosts the
3D topological Dirac state. The Dirac state has been described by the band
structure from first-principles calculations. The superconductivity was
characterized by the measurements of resistivity. The 2D nature of
superconductivity was verified by the angle dependent MR, supporting that the
superconductivity in BaMg2Bi2 occurs on the surface rather than in the bulk.

## Methods

###  0.1 Crystal growth and transport measurements

Large single crystals of BaMg2Bi2 were grown by flux method. The starting
elements, Ba, Mg and Bi, were mixed and put into the alumina crucible with a
molar ratio of Ba : Mg : Bi=1 : 5 : 9. The crucible was then sealed into an
evacuated quartz tube and heated to 900∘C. After cooling to 650∘C in 300 h,
the excess flux was removed with a centrifuge. The atomic composition of
BaMg2Bi2 has been checked by energy dispersive x-ray spectroscopy. The
measurements of transport properties were performed on a Quantum Design
physical property measurement system (PPMS). The obtained crystals are
sensitive to both air and GE varnish. The samples can survive in air for only
several minutes. The electrodes used in the measurements were made in a glove
box and then covered with a layer of N-grease.

###  0.2 Electronic structure calculations

The first-principles electronic structure calculations were performed with the
projector augmented wave (PAW) method[29, 30] as implemented in the Vienna _ab
initio_ simulation package (VASP)[31, 32, 33]. The generalized gradient
approximation (GGA) of Perdew-Burke-Ernerh (PBE) type[34] was adopted for the
exchange-correlation potential. The kinetic energy cutoff of the plane wave
basis was set to be 400 eV. For the Brillouin zone sampling, a
15$\times$15$\times$9 _k_ -point mesh was employed. The Gaussian smearing with
a width of 0.02 eV was used around the Fermi surface. In structural
optimization, both cell parameters and internal atomic positions were allowed
to relax until all forces on atoms were smaller than 0.01 eV/Å. The relaxed
structure parameters are well agreed with experimental results. The spin-
orbital-coupling (SOC) effect was included in the calculations of the
electronic properties.

## Data availability

The authors declare that the data supporting the findings of this study are
available within the article and its Supplementary Information. Extra data are
available from the corresponding authors upon reasonable request.

## References

## References

  * [1] Saito, Y., Nojima, T. & Iwasa, Y.  Highly crystalline 2D superconductors.  _Nat. Rev. Mater._ 2, 16094 (2017).
  * [2] Reyren, N. _et al._ Superconducting interfaces between insulating oxides.  _Science_ 317, 1196–1199 (2007).
  * [3] Wang, Q.-Y. _et al._ Interface-induced high-temperature superconductivity in single unit-cell FeSe films on SrTiO3.  _Chin. Phys. Lett._ 29, 037402 (2012).
  * [4] Ge, J.-F. _et al._ Superconductivity above 100 K in single-layer FeSe films on doped SrTiO3.  _Nat. Mater._ 14, 285–289 (2015).
  * [5] Liu, C. _et al._ Two-dimensional superconductivity and anisotropic transport at KTaO3 (111) interfaces.  _Science_ 371, 716–721 (2021).
  * [6] Liao, M. _et al._ Superconductivity in few-layer stanene.  _Nat. Phys._ 14, 344–348 (2018).
  * [7] Cui, J. _et al._ Transport evidence of asymmetric spin–orbit coupling in few-layer superconducting 1Td-MoTe2.  _Nat. Commun._ 10, 2044 (2019).
  * [8] Li, D. _et al._ Superconductivity in an infinite-layer nickelate.  _Nature_ 572, 624–627 (2019).
  * [9] Shen, D. _et al._ Two-dimensional superconductivity and magnetotransport from topological surface states in AuSn4 semimetal.  _Commun. Mater._ 1, 56 (2020).
  * [10] Song, J. _et al._ Coexistence of surface superconducting and three-dimensional topological Dirac states in semimetal KZnBi.  _Phys. Rev. X_ 11, 021065 (2021).
  * [11] Huang, C. _et al._ Proximity-induced surface superconductivity in Dirac semimetal Cd3As2.  _Nat. Commun._ 10, 2217 (2019).
  * [12] Xing, Y. _et al._ Surface superconductivity in the type II Weyl semimetal TaIrTe4.  _Natl. Sci. Rev._ 7, 579–587 (2020).
  * [13] Qin, S. _et al._ Quasi-1D topological nodal vortex line phase in doped superconducting 3D Dirac semimetals.  _Phys. Rev. Lett._ 123, 027003 (2019).
  * [14] Kobayashi, S. & Sato, M.  Topological superconductivity in Dirac semimetals.  _Phys. Rev. Lett._ 115, 187001 (2015).
  * [15] Qi, X.-L. & Zhang, S.-C.  Topological insulators and superconductors.  _Rev. Mod. Phys._ 83, 1057 (2011).
  * [16] Sato, M. & Ando, Y.  Topological superconductors: a review.  _Rep. Prog. Phys._ 80, 076501 (2017).
  * [17] Zhang, P. _et al._ Observation of topological superconductivity on the surface of an iron-based superconductor.  _Science_ 360, 182–186 (2018).
  * [18] Liu, Q. _et al._ Robust and clean Majorana zero mode in the vortex core of high-temperature superconductor (Li0.84Fe0.16)OHFeSe.  _Phys. Rev. X_ 8, 041056 (2018).
  * [19] Fu, L. & Kane, C. L.  Superconducting proximity effect and Majorana fermions at the surface of a topological insulator.  _Phys. Rev. Lett._ 100, 096407 (2008).
  * [20] Wang, M.-X. _et al._ The coexistence of superconductivity and topological order in the Bi2Se3 thin films.  _Science_ 336, 52–55 (2012).
  * [21] Xu, J.-P. _et al._ Experimental detection of a Majorana mode in the core of a magnetic vortex inside a topological insulator-superconductor Bi2Te3/NbSe2 heterostructure.  _Phys. Rev. Lett._ 114, 017001 (2015).
  * [22] Xu, J.-P. _et al._ Artificial topological superconductor by the proximity effect.  _Phys. Rev. Lett._ 112, 217001 (2014).
  * [23] Zhu, W. _et al._ Interfacial superconductivity on the topological semimetal tungsten carbide induced by metal deposition.  _Adv. Mater._ 32, 1907970 (2020).
  * [24] Zhu, W. _et al._ Superconducting interfaces between Weyl semimetal and normal metal.  _Adv. Quantum Technol._ 3, 2000020 (2020).
  * [25] Wang, H. _et al._ Discovery of tip induced unconventional superconductivity on Weyl semimetal.  _Sci. Bull._ 62, 425–430 (2017).
  * [26] Wang, H. _et al._ Observation of superconductivity induced by a point contact on 3D Dirac semimetal Cd3As2 crystals.  _Nat. Mater._ 15, 38–42 (2016).
  * [27] Aggarwal, L. _et al._ Unconventional superconductivity at mesoscopic point contacts on the 3D Dirac semimetal Cd3As2.  _Nat. Mater._ 15, 32–37 (2016).
  * [28] Takane, D. _et al._ Dirac semimetal phase and switching of band inversion in XMg2Bi2 (X= Ba and Sr).  _Sci Rep_ 11, 21937 (2021).
  * [29] Blöchl, P. E.  Projector augmented-wave method.  _Phys. Rev. B_ 50, 17953–17979 (1994).
  * [30] Kresse, G. & Joubert, D.  From ultrasoft pseudopotentials to the projector augmented-wave method.  _Phys. Rev. B_ 59, 1758–1775 (1999).
  * [31] Kresse, G. & Hafner, J.  Ab initio molecular dynamics for liquid metals.  _Phys. Rev. B_ 47, 558–561 (1993).
  * [32] Kresse, G. & Furthmüller, J.  Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set.  _Comp. Mater. Sci._ 6, 15–50 (1996).
  * [33] Kresse, G. & Furthmüller, J.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  _Phys. Rev. B_ 54, 11169–11186 (1996).
  * [34] Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).

This work is supported by the National Natural Science Foundation of China
(Nos.12104011, 11874336, 12074425, 11874422, 11904003, 12104010, 52102333,
12004003), the Natural Science Foundation of Anhui Province (Nos.2108085QA22,
1908085MA09, 2108085MA16), the Joint Funds of the National Natural Science
Foundation of China (No. U1832209) and the National Key R$\&$D Program of
China (No. 2019YFA0308602).

Y.-Y.W. coordinated the project and designed the experiments. Q.L. and Y.-Y.W.
synthesized the single crystals of BaMg2Bi2. P.-J.G. performed _ab initio_
calculations. Y.-Y.W. performed the transport measurements with the assistance
of Q.L. and X.-Y.Y.. Q.L., Y.-Y.W. and P.-J.G. plotted the figures and
analysed the experimental data. Y.-Y.W., X.-F.S. and T.-L.X. wrote the paper.
All authors discussed the results and commented on the manuscript.

accompanies this paper.

The authors declare no competing interests. The data that support the findings
of this study are available from the corresponding authors T.-L.X.
(tlxia@ruc.edu.cn), X.-F.S<EMAIL_ADDRESS>and Y.-Y.W<EMAIL_ADDRESS>upon reasonable request.

Figure 1: Crystal structure and electronic structure of BaMg2Bi2. a, Schematic
crystal structure of BaMg2Bi2. b,View of the structure from the direction of
_c_ -axis. c, The calculated band structure of BaMg2Bi2 along high symmetry
lines with the SOC effect included. d, The enlarged view of band structure
around $\Gamma$ point.

Figure 2: Superconductivity of the BaMg2Bi2 single crystals. a, Temperature
dependence of zero field resistivity $\rho$(_T_) in the range of 1.8 K to 300
K. Inset: the enlarged part from 1.8 K to 6 K. b, The $\rho$(_T_) curves under
different magnetic fields at low temperature. The applied current is 100
$\mu$A. c, The $\rho$(_T_) curves under various currents at low temperature
and zero magnetic field.

Figure 3: The critical field and critical current of the superconducting
state. a, Temperature dependence of critical field extracted from Fig. 2b. The
red line is a fit to the data using the empirical formula. b, Temperature
dependence of critical current obtained from Fig. 2c.

Figure 4: Angle dependent magnetoresistance in the superconducting and
nonsuperconducting states. a, d, Schematic of two configurations of the
measurements. The current is along the hexagonal edge and the magnetic field
rotates in the _yz_ -plane (a) and _xz_ -plane (d), respectively. b, e, Polar
plots of the angle dependent resistivity corresponding to the configurations
in a and d. The temperature is fixed at 1.8 K and the magnetic field varies
from 3 T to 14 T. c, f, Polar plots of the angle dependent resistivity
corresponding to the configurations in a and d. The magnetic field is fixed at
5 T and the temperature changes from 1.8 K to 12 K.


__

# Two-Phase Dynamics of DNA Supercoiling based on DNA Polymer Physics

Biao Wan  Beijing Computational Science Research Center, Beijing 100193,China
Jin Yu<EMAIL_ADDRESS>Department of Physics and Astronomy, Department of
Chemistry, NSF-Simons Center for Multiscale Cell Fate Research, University of
California, Irvine, CA 92697, USA

###### Abstract

DNA supercoils are generated in genome regulation processes such as
transcription and replication, and provide mechanical feedback to such
processes. Under tension, DNA supercoil can present a coexistence state of
plectonemic (P) and stretched (S) phases. Experiments have revealed the
dynamic behaviors of plectoneme, e.g. diffusion, nucleation and hopping. To
represent these dynamics with computational changes, we demonstrated first the
fast dynamics on the DNA to reach torque equilibrium within the P and S
phases, and then identified the two-phase boundaries as collective slow
variables to describe the essential dynamics. According to the time scale
separation demonstrated here, we accordingly developed a two-phase model on
the dynamics of DNA supercoiling, which can capture physiologically relevant
events across time scales of several orders of magnitudes. In this model, we
systematically characterized the slow dynamics between the two phases, and
compared the numerical results with that from the DNA polymer physics-based
worm-like chain model. The supercoiling dynamics, including the nucleation,
diffusion, and hopping of plectoneme, have been well represented and
reproduced, using the two-phase dynamic model, at trivial computational costs.
Our current developments, therefore, can be implemented to explore multi-scale
physical mechanisms of the DNA supercoiling-dependent physiological processes.

## INTRODUCTION

Supercoiling is ubiquitous in cellular DNA, which results from the topology of
double-helical DNA Bates  _et al._ (2005); Fuller (1971); Koster  _et al._
(2010); Lavelle (2014). For eukaryotic DNA, the large structural domains
($\sim 1$ M$b$) emerge in the genome as the topologically associating domains
(TADs)Dixon  _et al._ (2012); Glinsky (2015); Stadhouders  _et al._ (2019),
which are further divided into smaller supercoiling domains($\sim
10^{5}bp$)Naughton  _et al._ (2013); Gilbert and Allan (2014); Corless and
Gilbert (2016). For the prokaryotic DNA, comparatively large structural
domains and smaller supercoiling domains also present (e.g., $E$. $coli$ with
an average domain size, $\sim 10^{4}$ $bp$) Postow  _et al._ (2004); Gilbert
and Allan (2014). Supercoils are generated in many crucial genetic processes
and in turn regulates corresponding processes, including transcription,
replication and genome condensationLiu and Wang (1987); Ma  _et al._ (2013);
Chong  _et al._ (2014); Bates  _et al._ (2005); Burnier  _et al._ (2008);
Schvartzman  _et al._ (2013); Hirano (2000); Strick  _et al._ (2004). During
transcription elongation, a RNA polymerase (RNAP) moves along the helical DNA
and generates twists, which subsequently create positive supercoiling ahead of
the RNAP or downstream along the transcription direction and negative
supercoiling behind or upstreamLiu and Wang (1987). The accumulation of the
torsional stress can slow down or even stall the transcription elongationMa
_et al._ (2013); Chong  _et al._ (2014). Recent experiments show long-distance
cooperative and antagonistic dynamics of multiple RNAPs via supercoilingKim
_et al._ (2019), which justify the role of supercoiling as a long-range
mediatorTeves and Henikoff (2014). A previous modeling study pointed out that
multiple RNAPs facilitate elongation by reducing collisions or traffic jam by
torquesTantale  _et al._ (2016); and another model suggested that multiple
RNAPs also induce strong supercoiling to bring significant transcription
fluctuationsJing  _et al._ (2018). In addition, supercoiling contributes to
the chromosome condensationHirano (2000); Strick  _et al._ (2004), during
which the negative supercoils accumulate in the protein-free region of the DNA
Kimura and Hirano (1997). Under tension, the accumulative supercoiling
presents twisted and interwound coils called plectonemes, which may act as a
driving force for compacting the chromatin fiberKimura and Hirano (1997);
Swedlow and Hirano (2003).

Supercoil can be partitioned into two parts, twist and writhe, which play
distinct roles in genomic processesSeol and Neuman (2016). Underwound DNA
(with negative twists) can melt the secondary structures of DNA that can
attract regulatory proteinsMa and Wang (2014a, b); Seol and Neuman (2016). In
contrast, overwound DNA (with positive twists) can hinder enzymatic activities
associated with opening of the DNA duplex, such as in transcription initiation
and elongation by RNAPMa  _et al._ (2013); Chong  _et al._ (2014). On the
other hand, writhe characterizes spatial crossings of DNA segmentsFuller
(1971, 1978), accordingly facilitating distant interactions on DNAParker and
Halford (1991); Liu  _et al._ (2001); Huang  _et al._ (2001); Polikanov  _et
al._ (2007). Writhe can be in a solenoidal or a plectonemic form, or a
combination of both. Under an externally stretching force, supercoiled DNA can
present a coexistence state of bothMarko (1997a); Strick  _et al._ (1999);
Ghatak and Mahadevan (2005); Marko (2007), in which the plectonemic coils form
the plectonemic (P) phase, while the solenoidal coils are almost straightly
stretched, forming the stretched (S) phaseMarko (2007); Forte  _et al._
(2019).

DNA supercoiling has been quantitatively studied with single-molecule
techniquesMarko and Siggia (1994); Bustamante  _et al._ (1994); Strick  _et
al._ (1996, 1998, 1999); Bustamante  _et al._ (2000); Strick  _et al._ (2000);
Koster  _et al._ (2005); van Loenhout  _et al._ (2012); Bustamante  _et al._
(2003). The extension-twist curves from experiments reveal a coupling between
supercoiling (linking number) and DNA extension under stretching forceMarko
and Siggia (1994); Bustamante  _et al._ (1994); Strick  _et al._ (1996). The
discontinuities on torque and extension versus linking number were also
observed in experimentForth  _et al._ (2008). Experiments with a magnetic
tweezers pulling a fluorescently labeled DNA have directly imaged the
coexistence phase, and remarkably, measured the time-dependent supercoiling
dynamics, e.g., plectoneme diffusion, nucleation and hoppingvan Loenhout  _et
al._ (2012). The propagation of plectoneme along DNA via diffusion presents
conformational rearrangement of DNA at about hundreds bps per second. The
hopping happens comparatively fast and a fastest one observed takes tens of
milliseconds for a displacement of thousands of bpsvan Loenhout  _et al._
(2012). The surviving time of the individual plectoneme from nucleation to
vanishing, i.e., the plectoneme life-time, spans from milliseconds to several
seconds.

In accompany with the experiments, theoretical studies of DNA supercoiling
have been developed based on the worm-like chain (WLC) modelMarko and Siggia
(1995a, b); Moroz and Nelson (1997); Marko (2007); Phillips  _et al._ (2009).
The statistical mechanics of the WLC model of DNA provides a description of
the coupling between supercoiling and DNA extensionMoroz and Nelson (1997);
Marko (1998), and reveals the pre-plectonemic loops formed by bending DNA,
i.e., buckling transitionDaniels and Sethna (2011), which give rise to the
discontinuities on torque and extension versus linking number curvesMarko and
Neukirch (2012); Emanuel  _et al._ (2013). Further, studies show that under a
stretching force, the coexistence of the S and P phases is maintained under a
coexistence torque that is a function of the stretching forceMarko (2007).

Nevertheless, the related analytic theories concern mostly the equilibrium and
static behaviors, without elucidating the time-dependent dynamic processes.
The experiments suggest that the dynamics of plectonemes essentially covers
multiple orders of timescalesvan Loenhout  _et al._ (2012). In this study, we
focus on building a model to describe these dynamic processes across the time
scales.

Apart from the analytic theories, computational simulations and numerical
methods have offered effective ways to probe supercoiling dynamics. Atomistic
molecular dynamics (MD) simulations offer the finest details of such dynamics.
For example, MD simulations of DNA mini-circles ($10^{2}$ $b$) reveal writhe
fluctuation and configuration diversityMitchell and Harris (2013); Mitchell
_et al._ (2011); Irobalieva  _et al._ (2015). To improve computational
efficiency, a coarse-grained model that treats nucleotides as beads with three
interaction sites, called oxDNA, was developedMatek  _et al._ (2015).
Furthermore, by unitizing MD and the coarse-grained simulations, the multi-
scale dynamics of supercoiling have been studiedSutthibutpong  _et al._
(2016); Hirsh  _et al._ (2013a). Nevertheless, computational costs of these
models are high and the simulation time scales are limited by microsecondsDans
_et al._ (2016). The models without considering DNA sequence structure are
capable of exceeding such a limit. An elasto-dynamic model (i.e., the
Kircholff rod), for example, has been developed, by which DNA is characterized
in terms of a continuum rod while omitting thermal fluctuationsGoyal  _et al._
(2005). The corresponding applications were carried out on supercoil
removalLillian  _et al._ (2011), and on compressed DNA inside bacteriophage
cavity to allow DNA ejectionHirsh  _et al._ (2013b). In particular, for a
representative polymer physics model, the numerical approach of the WLC model,
i.e., the discrete worm-like chain (dWLC) method, has been widely implemented
in simulating DNA with the thermal fluctuations incorporatedVologodskii
(2006). For example, the Monte Carlo (MC) simulations of the WLC model have
been used in describing DNA thermodynamics, conformation and site
juxtapositionVologodskii  _et al._ (1992); Klenin  _et al._ (1991); Polikanov
_et al._ (2007). The Brownian dynamics (BD) simulations, on the other hand,
have been applied more widely in studying supercoilingVologodskii (2006);
Klenin  _et al._ (1998). Such type of studies successfully show the plectoneme
diffusionBell  _et al._ (2012), supercoiling conformations Vologodskii  _et
al._ (1992); Babamohammadi and Lillian (2020), buckling transitionDaniels and
Sethna (2011); Walker  _et al._ (2018); Ott  _et al._ (2020), modeling a DNA
wrapped around a model histoneWocjan  _et al._ (2009), and supercoil removal
by nickingIvenso and Lillian (2016). Besides, other semi-flexible polymer
models with Lennard-Jones potential have also been utilized to investigate the
mechanism for supercoiling by rotating the endsWada and Netz (2009), or using
interwoven and braided polymers to produce plectonemesForte  _et al._ (2019).

The computational costs of the above models are still quite significant for
studying the plectoneme dynamics. For example, by employing the dWLC method,
simulating the plectoneme dynamics over a second requires weeks of CPU
timeIvenso and Lillian (2016); Ott  _et al._ (2020). In this study, we
developed a two-phase (S and P phases) dynamic model of DNA supercoiling at
trivial computational cost based on the WLC polymer physics model of DNA. In
solvent, supercoiled DNA experiences frictional forces on the two phases and
the energy gradient that drives the transformation between the S and P phases.
Thus the first thing is to identify the fast dynamics comparing to the phase-
transformation and possibly integratae them out to reduce the degrees of
freedomFasnacht  _et al._ (2004); Zhang and Ma (2010). Here we demonstrate
first the dynamics and timescales of the torque transport within S and P
phases that define the fast dynamics. Then we choose interphase boundaries as
the slow dynamic observables and obtain the associated energetics.
Subsequently, we derive the Langevin dynamics of the slow observables. As a
result, this implies that DNA supercoil can fast propagate as twist (torque)
and slowly propagate as writhe (or plectonemes) along DNA. For calibration and
consistency check, we compare the numerical results of the two-phase dynamic
model with those from the dWLC method. For applications, we reproduce the
supercoil dynamics on plectoneme diffusion, nucleation and vanishing as being
measured experimentally, as well as the force/ionic strength-dependent
dynamics. Finally, plectoneme hopping and the associated energetics
transformations between the two phases are monitored.

## RESULTS

Here, we study dynamics of the DNA supercoiling and build a two-phase dynamic
model with a focus on DNA of $10^{3}$ to $10^{4}$ $b$. Such sizes of DNAs are
widely used in the single-molecule experimentsCrut  _et al._ (2007); van
Loenhout  _et al._ (2012); Ma  _et al._ (2013). DNAs shorter than several
hundreds bps hardly bend to form a stable plectonmic helix; for DNAs longer
than $10^{4}$ $b$, however, the supercoiling domain may emerge and take
overDixon  _et al._ (2012); Glinsky (2015).

### Fast torsional equilibrium on the stretched (S) and plectonemic (P) phases

In this section, we describe the dynamics of torque transport reaching to
equilibrium within the S and P phases, and show that such dynamics are fast
comparing to the transformation between the two phases.

#### Torque transport on stretched DNA or S phase

We study the dynamics of torque propagation on supercoiled DNA based on the
WLC model numerically, i.e., via an implementation of the dWLC method (see
Text S1 in the Supporting Material).

In the WLC model, harmonic potentials are used in representing the
elasticities of DNABates  _et al._ (2005); Marko and Siggia (1995b); Klenin
_et al._ (1998). The twist energy density is defined as

$\mathcal{E}_{twist}=\frac{1}{2}k_{B}Tl_{t}(\frac{\Delta\theta}{\Delta
u})^{2}$ (A1)

where $k_{B}$ is the Boltzmann constant, $T$ is the room temperature, $l_{t}$
is the torsional persistence length, i.e., $l_{t}=75nm$, $\Delta\theta$ is the
twist on the segment $\Delta u$. We introduce an internal torque $T_{q}^{in}$
as a measure of torsional stress on DNA, $T_{q}^{in}\equiv
k_{B}Tl_{t}\frac{\Delta\theta}{\Delta u}$ (Torque is proportional to twist).
Based on the elastic properties of DNA, the torque transport along DNA can be
written as (see Text S2),

$\frac{\partial T_{q}^{in}}{\partial
t}=\frac{l_{t}k_{B}T}{\zeta_{R}}\frac{\partial^{2}T_{q}^{in}}{\partial u^{2}}$
(A2)

where drag coefficient for a rotational cylinder $\zeta_{R}=4\pi R^{2}\eta$,
and the radius of DNA $R=2nm$ and the viscosity $\eta=0.001kg/m/s)$. A
solution to Eq A2 can be seen in Text S2. The characteristic time for torque
equilibrium time can be reached at
$2.5\frac{\zeta_{R}}{l_{t}k_{B}T}L^{2}$(about 0.1$ms$ for $3000bp$) (Text S2),
independent of stretching force, torque and linking number. For example, a
twist pulse starting at one end of the stretched DNA or S phase will spread
over $10^{3}$ $bp$ within $10^{-2}ms$ and $10^{4}$ $bp$ within $1ms$. Compared
with the dynamics of plectoneme of $\sim 10^{3}$ $bp$ that happens at the time
scale from $ms$ to $s$, e.g. diffusion and hopping, the torque transport
within the S phase is fast.

Figure 1: Measuring the torque transport on the stretched (S) phase and the
plectonemic (P) phase of the supercoiled DNA based on the dWLC, respectively.
(a)The schematic shows the torque transport along the S phase following an
external torque $T_{q}^{ex}$ promptly applied from 0 $pNnm$ to 9 $pNnm$ at
$t=0$. (b) The internal torque of DNA at length of $1500bp$ and $3000bp$,
$\bar{T}_{q}^{in}$ rises to equilibrium (9 $pNnm$) within 0.1 $ms$ to follow
$T_{q}^{ex}$. (c) The excess linking number density $\sigma_{S}$ reaches to
equilibrium $\sigma_{S}^{eq}=0.17$ within 0.1$ms$. (d) The schematic shows
that a pure plectonemic DNA follows an promptly applied torque,$T_{q}^{ex}$
from an initial $8pNnm$ to maintains the plectonemic helix ($t<0$) to $10pNnm$
at $t=0$. The internal torque along DNA consequently propagates on plectoneme.
(e) The mean torque of DNA at length of $1500bp$ and $3000bp$,
$\bar{T}_{q}^{in}$ reaches to equilibrium (10 $pNnm$) within $0.3ms$. (f) The
excess linking number density $\sigma_{P}$ reaches to equilibrium value
$\sigma_{P}^{eq}=0.091$ within 0.3 $ms$.

Fig.1 shows the scheme for simulating the torque transport on a stretched DNA
using the dWLC method. Under a constant stretching force $f=0.4$ $pN$, DNA is
twisted by an external torque $T_{q}^{ex}=9pNnm$ starting at $t=0$
($T_{q}^{ex}=0$ at $t<0$). After that, the torque is kept constant and the
torsional response of DNA is monitored.

The measurements of the mean internal torque of DNA at lengths of 1500 $bp$
and 3000 $bp$ are shown in Fig. 1. After being initially twisted, the mean
internal torque, $\bar{T}_{q}^{in}$, reaches to equilibrium (9 $pNnm$) at
about $0.03ms$ for 1500 $bp$ and about $0.1$ $ms$ for 3000 $bp$, consistent
with the torque equilibrium time $2.5\frac{\zeta_{R}}{l_{t}k_{B}T}L^{2}$.

Supercoiling can be quantified by the excess linking number $\Delta Lk$ or the
linking number density $\sigma$ (see Text S3). The free energy density of a
stretched DNA is,
$\mathcal{F}_{S}=\frac{1}{2}c_{S}(\sigma_{S}^{eq})^{2}+\mathcal{F}_{0}$(Text
S4), including the contributions from torsion and force stretchingMarko (1998,
2007), where $c_{S}$ stands for the torsional stiffness, which is a function
of $f$, $l_{t}$ and $l_{b}$ ($l_{b}=50nm$, the DNA persistence length),
$\sigma_{S}^{eq}$ is the equilibrium linking number density of the stretched
phase, and$\mathcal{F}_{0}$ is the free energy density without torsion,
representing the stretching of the WLCMarko and Siggia (1995a). At
equilibrium, the torque is defined as
$\frac{1}{\Omega_{0}}\frac{\partial\mathcal{F}_{S}}{\partial\sigma_{S}^{eq}}$Marko
(2007), where the constant $\Omega_{0}=1.85nm^{-1}$ is the rotation angle of
the DNA backbone per unit lengthBates  _et al._ (2005); Marko (2007). Thus, we
have

$\sigma_{S}^{eq}=\frac{\Omega_{0}\bar{T}_{q}^{in}}{c_{S}}$ (A3)

We accordingly measure the excess linking number density of the stretched DNA
as shown in Fig. 1. $\sigma_{S}$ reaches to the equilibrium value at 0.03$ms$
for 1500 $bp$ and 0.1 $ms$ for 3000 $bp$.

#### Torque transport on plectoneme

A plectoneme can be considered as a pair of rods anti-parallel interwound
around each other under an external torqueRybenkov  _et al._ (1997); Marko
(2007). Thus, an empirical form of the energy of the plectonemeKlenin  _et
al._ (1991); Rybenkov  _et al._ (1997); Marko (2007) can be utilized,

$\mathcal{F}_{P}=\frac{1}{2}c_{P}\sigma_{P}^{2},$ (A4)

where $c_{P}\equiv l_{p}k_{B}T\Omega_{0}^{2}$ is the torsional stiffness of
plectoneme, and $l_{p}$ is the torsional persistence length of the plectonemic
helix, depending on the ionic strength and often taking value of 21 to 27 $nm$
Rybenkov  _et al._ (1997), and $\sigma_{P}$ is the excess linking number
density of the plectonemes. The internal toque is defined similarly as that of
the stretched DNA, i.e., $T_{q}^{in}=k_{B}Tl_{p}\sigma_{P}$. Nevertheless,
torque transport on plectoneme can not be described as the twist angle
$\theta$ propagation on DNA. That is because that the change of twist of the
plectomeical helix also influences the writhe, and in turn the writhe change
feeds back to the twist propagation. We introduce a plectonemic angle $\Theta$
(see Fig. S1), similar to $\theta$, as
$T_{q}^{in}=k_{B}Tl_{p}\frac{\Delta\Theta}{\Delta u}$. Thus we consider the
torque transport on plectoneme as the propagation of the plectonemical angle,
i.e.,

$\frac{\partial T_{q}^{in}}{\partial
t}=\frac{l_{p}k_{B}T}{\zeta_{P}}\frac{\partial^{2}T_{q}^{in}}{\partial
u^{2}},$ (A5)

where $\zeta_{P}\approx 2\pi R_{plect}^{2}\eta$ is the drag coefficient due to
motions of the two parallel rods and $R_{plect}$ the radius of
plectonemesNeukirch and Marko (2011); van Loenhout  _et al._ (2012), about 2.5
nm. The time for torque equilibrium on plectoneme with contour length $L$ can
be calculated as $2.5\frac{\zeta_{P}}{l_{p}k_{B}T}\frac{L^{2}}{4}$, about
$2.6$ times comparing to that on the stretched DNA. i.e, about 0.26 $ms$ for
torque equilibrium for 3000 $bp$.

Fig. 1 illustrates the setup of torque transport on a pure plectoneme, which
is originally maintained by a pre-existing external torque
$T_{q}^{ex}=8pNnm$($t<0$). Then the external torque jumps to
$T_{q}^{ex}=10pNnm$ at $t\geqslant 0$, and the ensuing torsional response of
the plectoneme is monitored.

The measurements of the internal torque of plectonemes at lengths of 1500 $bp$
and 3000 $bp$ are shown in Fig. 1. Following the torque jumping at $t=0$, the
mean internal torques, $\bar{T}_{q}^{in}$ rises from $8pNnm$ to $10pNnm$ at
about 0.1 and 0.3 $ms$ for 1500 $bp$ and 3000$bp$, respectively, as estimated
above.

In equilibrium, Eq. A4 becomes
$\mathcal{F}_{P}=\frac{1}{2}c_{P}(\sigma_{P}^{eq})^{2}$. The equilibrium
torque is defined as
$\frac{1}{\Omega_{0}}\frac{\partial\mathcal{F}_{P}}{\partial\sigma_{P}^{eq}}$.
Thus

$\sigma_{P}^{eq}=\frac{\Omega_{0}\bar{T}_{q}^{in}}{c_{P}}$ (A6)

We also numerically measure the excess linking number density of the
plectoneme as shown in Fig. 1. As Eq A6 suggests,
$\sigma_{P}^{eq}=\frac{T_{q}^{ex}\Omega_{0}}{c_{P}}=0.076$ at $t<0$, where
$T_{q}^{ex}=8pNnm$. Then $\sigma_{P}$ follows the external torque jumping.
Subsequently, $\sigma_{P}$ reaches a new equilibrium defined by
$\frac{T_{q}^{ex}\Omega_{0}}{c_{P}}=0.091$ at about $0.1$ ms for 1500 $bp$ and
0.3 ms for 3000 $bp$.

The torque equilibrium (Eqs A3 and A6) can be applied to the coexistence state
of S and P phases. The interphase torque equilibrium means the equilibrium at
the phase-boundaries, i.e., $c_{P}\sigma_{P}^{eq}=c_{S}\sigma_{S}^{eq}$. Next,
we show that the phase boundaries can be chosen as the slow observables to
quantify the plectoneme dynamics.

### Two-phase dynamics of DNA supercoiling

We have demonstrated the comparatively fast dynamics of torque to reach
equilibrium within the S and P phases, respectively. The fast dynamics can
then be averaged out within the two phases. We show below that the energy
gradients at the phase-boundaries drive the transformation between the two
phases and thus establish a two-phase dynamic model for describing the DNA
supercoiling.

Multiple plectonemes have been observed in the extended DNA supercoil in
experimentvan Loenhout  _et al._ (2012). The dynamics of the plectonemes can
be described by the boundaries between the S and P phases. We therefore label
the plectonemes from the fixed end to the stretching end with
$\alpha=1,2,...$. Without loss of generality, we only focus on the $\alpha$-th
plectoneme as shown in Fig. 2.
${\mathbf{X}}\equiv(X_{\alpha}^{l},X_{\alpha}^{r})$ are the boundaries of
plectonemes, where the superscripts $l$ and $r$ denote the left and right
boundaries of the $\alpha$-th plectoneme, respectively. Along the axis,
$X_{\alpha}^{l}<u\leqslant X_{\alpha}^{r}$ locates the $\alpha$-th plectoneme,
while $X_{\alpha-1}^{r}<u<X_{\alpha}^{l}$ specifies the $\alpha$-th section of
the S phase. Thus the motion of plectonemes must be accompanied by that of the
stretched phase. The formation of plectonemic coils can make DNA compact by
bringing distant-site together, i.e., $X_{\alpha}^{l}$ and $X_{\alpha}^{r}$
are in physical contact. Plectoneme propagation corresponds to the movement of
$X_{\alpha}^{l}$ and $X_{\alpha}^{r}$ in the same direction. For example, in
Fig. 2 if the $\alpha$-th plectoneme moves left, the S phase on the left-hand
side of $X_{\alpha}^{l}$ is transformed into the P phase while the P phase on
the left-hand side of $X_{\alpha}^{r}$ is transformed into the S phase. This
process is accompanied with the internal slithering of parallel segments of
the plectonemic helix(Fig. 2). Indeed, the velocity of the relative slithering
is equal to the velocity of plectoneme motion (see Methods).

Figure 2:  Two-phase dynamic model of DNA supercoiling. (a) The schematic of
the $\alpha$-th plectoneme among multiple plectonemes, showing the torque
equilibrium between the plectoneme and stretched phases. By mapping the P and
S phases into a 1-dimensional arc-length parameter space, one introduces
$u=X_{\alpha}^{l}$, the left boundary of the $\alpha$-th plectoneme, and
$u=X_{\alpha}^{r}$, the right boundary of the $\alpha$-th plectoneme. (b) The
propagation of the $\alpha$-th plectoneme. The slithering of the segments of
the plectonemic helix facilitates the recombination of distal sites on DNA
(marked by red and black marks). (c)Simulation trajectories based on the WLC
model and the two-phase dynamics. The internal torque of a DNA of 4500 $bp$ at
$f=0.5pN$ fluctuates around the coexistence torque, $T_{co}$. The histograms
in the inset show the variance of the internal torque obtained from the dWLC
simulations (red bars, 50 $ms$ in total) and that from the simulations using
the two-phase dynamic model (green bars, 1 $s$ in total).

Such observables always evolve slowly comparing to the torque transport,
independently of the contour-length (Text S5). A collective observable
$\mathbf{X}$ is associated with a certain amount of free energy, with its
gradient along $\mathbf{X}$ as a driving force to reach the free energy
minimaWeinan (2011). Based on this idea we derive below the free energy and
the dynamic equations associated with the slow variable $\mathbf{X}$.

####  Free energy associated with the slow observables

As the torque equilibrium is reached, i.e.,
$\Omega_{0}\overline{T}_{q}^{in}=c_{P}\sigma_{P}^{eq}$. The torque equilibrium
at the boundaries means $c_{P}\sigma_{P}^{eq}=c_{S}\sigma_{S}^{eq}$.
Accordingly, every plectenome carries the equal $\sigma_{P}^{eq}$. The free
energy density of each plectoneme is therefore represented as
$\mathcal{F}_{P}(\mathbf{X})=\frac{1}{2}c_{P}(\sigma_{P}^{eq})^{2}$Klenin  _et
al._ (1991); Rybenkov  _et al._ (1997).

The translational displacement of the $\alpha$-th section of the S phase is
subjected to the viscous drag along it. At $u<X_{\alpha}^{r}$ and
$>X_{\alpha}^{l}$ , the tension imposed on DNA is $f_{u}$, and at the
stretching end, $f|_{u=L}=f$. As the torque equilibrium is reached, the local
free energy density (or force) is represented by
$\mathcal{F}_{S}(u,\mathbf{X})=\frac{1}{2}c_{S}(\sigma_{S}^{eq})^{2}+\mathcal{F}_{0}$Marko
(1998, 2007).

The free energy associated with ${\mathbf{X}}$ can be written as the
contributions from the S and P phases

$\Phi_{0}(\mathbf{X})=\sum_{\alpha=0}\int_{X_{\alpha}^{r}}^{X_{\alpha+1}^{l}}\mathcal{F}_{S}(u,\mathbf{X})du+\int_{X_{\alpha+1}^{l}}^{X_{\alpha+1}^{r}}\mathcal{F}_{P}(\mathbf{X})du.$
(B1)

#### Langevin dynamics

The free energy gradient along $\mathbf{X}$ serves as a driving force, which
drives the transformation between S and P phases. The variable $\mathbf{X}$
thus follows the Langevin equations
$\gamma_{1}\mathbf{\dot{X}}=-\nabla_{\mathbf{X}}\Phi_{0}(\mathbf{X})+\sqrt{\frac{2k_{B}T\gamma_{1}}{dt}}\mathbf{\dot{W}}$,
where $\gamma_{1}$ denotes the drag coefficent for the growth or shrinkage of
plectonemes, and $\dot{W}$ is the white noise. Specifically, the left/right
boundaries of the $\alpha$-plectoneme follow

$\dot{X}_{\alpha}^{l/r}=-\frac{1}{\gamma_{\alpha}}\nabla_{X_{\alpha}^{l/r}}\Phi_{0}(\mathbf{X})+\sqrt{\frac{2k_{B}T}{\gamma_{\alpha}}}\dot{W},$
(B2)

where $\gamma_{\alpha}$ is drag coefficient for the growth or shrinkage of the
$\alpha$-th plectoneme (see Methods).

It should be noted that Eq B2 should be completed with two constraints: the
first one is the topology constraint ($\Delta Lk=\Delta Lk_{S}+\Delta
Lk_{P}.$), and the second one is in-extensibility of DNA model
($L=L_{S}+L_{P}$)(see Methods).

The equilibrium point, i.e., $\nabla_{X_{\alpha}^{l/r}}\Phi_{0}(\mathbf{X})=0$
leads to the coexistence torque
$T_{co}=\frac{1}{\Omega_{0}}(\frac{2c_{P}(f-\sqrt{\frac{k_{B}Tf}{l_{b}}})}{1-c_{P}/c_{S}})^{\frac{1}{2}}$Marko
(2007). Actually, the dynamics of the slow variable $\mathbf{X}$ accompanies
with a time-dependent torque,
$T_{q}^{in}(t)\equiv\frac{c_{P}\sigma_{P}^{eq}(t)}{\Omega_{o}}$ (Fig. 2). The
internal torque of a DNA at length of 4500 $bp$ under $f=0.5pN$ fluctuates
around the coexistence torque, $T_{co}$. The two-phase dynamic model can be
regarded as a smoothed version of the dWLC method. Although the fluctuations
faster than the torque equilibrium have been averaged out, the mean and the
variance of the torque are still consistent with that obtained from the dWLC
method(inset in Fig. 2).

Supercoil with a given linking number is analogous to the classical gas (van
der Waals) confined inside a fixed volume (NVT ensemble)Mattis and Swendsen
(2008). A DNA twisted under a constant torque $T_{q}^{ex}$, however, is
analogous to the classical gas under a fixed pressure (isobaric ensemble). Its
free energy densities can be obtained through the Legendre’s transformation
$\mathcal{F}_{S/P}(T_{q}^{ex})=\mathcal{F}_{S/P}-\Omega_{0}T_{q}^{ex}\sigma^{eq}$.
The free energy densities of the plectonemes and the stretched phase are
therefore
$\mathcal{F}_{P}=-\frac{1}{2}\frac{(\Omega_{0}T_{q}^{ex})^{2}}{c_{P}}$ and
$\mathcal{F}_{S}=-\frac{1}{2}\frac{(\Omega_{0}T_{q}^{ex})^{2}}{c_{S}}-f_{u}+\sqrt{\frac{k_{B}Tf_{u}}{l_{b}}}$,
respectively. Inserting them into Eq B1 and Eq B2 we obtained the dynamic
model of extended DNA supercoiling under constant torque.

Generally speaking, the stretching force, $f$ in Eq B1 or B2 is not
necessarily constant. The force can be time-dependent as long as it dose not
break the interphase torque equilibrium. For example, one can impose an
instantaneous force. Similarly, $\Delta Lk$ can be time-dependent as well if
it does not break the interphase torque equilibrium.

### Comparing with the WLC model

In this section, by employing the dWLC method, we compare and calibrate the
two-phase dynamic model via two examples. One is the supercoiling accumulation
at a constantly twisting rate. The other one is supercoiling accumulation
under a constant torque. In the first example, the system reaches to a two
phase-coexistent state at equilibrium. We show that the two-phase dynamic
model at trivial computational costs provides consistent results with the dWLC
method, which demands a significant amount of computational time. In the
second case, plectonemes finally dominate the supercoiling state under an
external torque larger than $T_{co}$. This non-equilibrium process helps
determine the slithering drag coefficient by calibration.

#### Supercoil accumulation and buckling transition at a constant rate

An extended DNA undergoes linking number accumulation when one end of the DNA
is rotated while the other end is kept fixed, as shown in Fig.3. The
stretching force $f=0.5pN$, the rotation angle $\Omega=\omega t$ and the
angular velocity $\omega=20\pi/s$, equivalent to $\Delta Lk=10$ per second
(i.e., comparable to polymerase enzyme unwinding and synthesizing at 100
nucleotide or nt per second)Thomen  _et al._ (2005); Jing  _et al._ (2018). At
first, the extended DNA is twisted, and then it undergoes a buckling
transition (discontinuities on torque and extension curves) followed by
plectoneme formationForth  _et al._ (2008).

We performed BD simulations of such a supercoil accumulation based on the dWLC
method. In Fig. 3, the relative extension $z/L$ almost remains unchanged
before 0.6 s while the torque gradually builds up, suggesting that the
extended DNA is twisted. The buckling at about 0.6 s marks the starting-point
of the phase-coexistence. The buckling transition indicates a competition
between pre-plectonemic loops and stretched phaseMarko and Neukirch (2012);
Mielke  _et al._ (2004). After 0.6 s, the torque drops and converges as
torsional stress on DNA is relieved, indicating the twist is transformed into
writhe. These discontinuities have previously been observed in single-molecule
experimentForth  _et al._ (2008). The coexistence torque $T_{q}^{co}$ remains
unchanged when more linking number is injected. Note that multiple plectonemes
can co-exist in the process. The representative snapshots are also shown in
Fig. 3. The unchanged torque curves or the linearly decreased extension after
buckling transition suggest that the process remains quasi-static.

It should be noted here that we obtained the curves in Fig. 3 from four BD
simulation trajectories using the dWLC method, and each of the trajectories
requires several weeks of CPU time. The calculation of the electrostatics
interaction between the charged points on segments of the dWLC model indeed
dominates the computational cost (see Text S6). According to our analyses, the
computational cost of the dWLC in the plectonemic phase is still significant
and grows linearly with the contour-length. For example, in Fig. 1, simulating
a stretched DNA of 3000 $bp$ in 10 $ms$ costs 3 hours while simulating a pure
plectoneme of the same contour length costs one day(i.e., the CPU time). The
two-phase dynamic model, however, enables us to reduce the computational cost
to trivial by dealing only with the essential variables with slow dynamics. In
Fig. 3, the torque and extension curves generated from the two-phase dynamic
model (using tens of seconds of computation time) well reproduce those in
Fig.3(b) from the dWLC method and also consistent with previous experimental
and computational studiesForth  _et al._ (2008); Ott  _et al._ (2020).

Figure 3: Buckling transition and measurements of the extensions of the DNA
during supercoil accumulation. (a)Injecting supercoil (i.e., increasing
linking number) at constant rate to the DNA of 3000 $bp$ within 1 $s$. The
rotation rate $\omega$ is 20 $\pi$ per second. (b)The relative extension $z/L$
(black curve) and torque (red curve) obtained from the BD simulations based on
the dWLC as the DNA is gradually twisted. (c)The torque and extension curves
from (b) are reproduced using the two-phase dynamic model. (d) Injecting
supercoil under a constant external torque, $T_{q}^{ex}=9pNnm$. (e)The
relative extension from the BD simulations based on the dWLC with various DNA
lengths. (f) The results from (e) are reproduced using the two-phase dynamic
model. The curves in (e) and (f) are the fitting functions $\sqrt{1+\lambda
t}$ (see Text S7).

#### Supercoil accumulation under a constant torque

We also performed simulations of supercoiling generation under a constant
torque. Fig.3 illustrates an extended supercoil under a constant torque,
$T_{q}^{ex}=9pNnm$ (comparable to 5$\sim$11 $pNnm$ torque exerted by
polymerasesMa  _et al._ (2013)), and a constant stretching force $f=0.4pN$.
Since the external torque is larger than $T_{co}$ (i.e., $6.5pNnm$ that
maintains the coexistence state under 0.4 $pN$), plectonemes are created and
the ultimate equilibrium state is dominated by plectonemes. Consequently,
plectonemic coils are building up until all stretched DNA is interwound to
plectonemical coils.

Using the dWLC method, we can capture the dynamics of the plectoneme growth.
Fig. 3 shows the measurements of the relative extensions $z/L$ of DNAs at 1500
$bp$ (black solid squares), 3000 $bp$, 6000 $bp$ and 9000 $bp$. Indeed, the
shrinking extensions indicate the growth of plectonemes. The fitting curves
are generated using $\frac{z}{L}-\frac{z_{0}}{L}\propto 1-\sqrt{1+\lambda t}$,
where $\lambda$ is a fitting parameter. In the two-phase dynamic model, all
but the slithering drag coefficients are given based on the hydrodynamics (see
Methods). The slithering drag $\mu_{slt}$ is to be determined in the two-phase
model by calibrating results with the dWLC method. We performed the
simulations using the two-phase dynamic model shown in Fig. 3. The slithering
drag $\mu_{slt}$ was then determined via fitting parameter $\lambda$ Text S7).
Consequently, we can then reproduce the plectoneme dynamics.

### Reproducing experimentally measured plectoneme dynamics

Previous experimental studies show that the plectonemes on the extended
supercoiled DNA can diffuse, vanish and re-emerge, and their dynamics depend
on the stretching force and ionic strengthvan Loenhout  _et al._ (2012).

Figure 4: Obtaining plectoneme dynamics using the two-phase dynamic model. (a)
The kymograph of the plectonemes of a 21-kbp (7$\mu m$) under an constant
stretching $f=0.8pN$ and with a constant linking number $\Delta Lk=70$. (b)The
diffusion coefficients of a single plectoneme under different forces obtained
from the simulation of the two-phase model (the dots). The green curve is
directly from the Einstein relation. (c) Large stretching force $f$ makes more
compact plectonemic helices (smaller radius of the plecotnemic helix). (d)
Increased salt concentration makes more compact plecotnemic helices. (e)
Buckling transition. A portion of stretched DNA is bent to form a loop during
nucleation, which is associated with an energy barrier. The reversed process
(red arrow, a plectoneme degenerates to a loop) to the nucleation (blue arrow)
is also illustrated. (f) Nucleation rate of the plectonemes as a function of
force obtained from the two-phase model simulations. The black, red and blue
curve and error bars are the results for ($l_{p}=25nm$), ($l_{p}=23nm$) and

($l_{p}=21nm$) corresponding to the low, regular and high salt concentrations,
respectively. (d) The mean lifetime of the individual plectoneme as a function
of force and salt concentration. Note that only a single plectoneme is
captured at $l_{p}=21nm$ and $f=3.2pN$ in simulations (nucleation and
annihilation, however, need multiple plectonemes). (h) A typical hopping event
of a displacement of plectoneme for nearly 4 $\mu m$ in a time interval about
$0.1s$, corresponding to 12 kbp. (i) Energetic conversion during hopping (
nucleation is, e.g, at 80 $ms$). The total energy changes of hopping (black
line) and the reversed (complementary to hopping) process (dashed line) are
also shown in the inset.

In particular, the experiments also captured plectoneme hopping, in which a
plectoneme vanishes at one position and a new one immediately emerges by
nucleation at another remote position. By employing the two-phase dynamic
model, we reproduced the plectoneme dynamics, including plectoneme diffusion,
nucleation and hopping. Notably, we captured the plectoneme hoppings with
significant energetic changes, which facilitate such rapid processes.

#### Reproducing plectoneme diffusion

The inter-phase boundaries in the two-phase dynamic model can be used to
directly monitor the plectoneme dynamics. Fig. 4 shows a kymograph, i.e., the
spatial vs time evolution process of the plectonemes of a 21-kbp (7$\mu m$)
DNA (extended under a force $f=0.8pN$ and with $\Delta Lk=70$), illustrating
the coexistence of multiple plectonemes, and the plectoneme nucleating and
vanishing.

The simulations of a plectoneme with $1/4$ of the DNA contour length under
different forces were performed. To estimate the diffusion constant, we
measured the mean square displacement (MSD) $\overline{\Delta x^{2}(t)}$ of
the plectoneme. As shown in Fig. 4, we obtained the diffusion constants via
$D\equiv\frac{\overline{\Delta x^{2}(t)}}{2t}$ under different stretching
forces. For consistency check, we also compared the results with the Einstein
relation, i.e., $D=\frac{k_{B}T}{\zeta_{diff}}$ where $\zeta_{diff}$ as a sum
of the viscous drags for the transverse displacement of the plectoneme (see
Methods). The diffusion constant of the plectoneme is about $10^{-1}\mu
m^{2}/s$.

#### Reproducing nucleation and Measuring plectoneme life time

To compare the dynamics of the plectonemes at different stretching forces and
ionic strengths, simulations of the P phase with $1/4$ of the DNA contour
length were performed in the conditions where the stretching forces were
$0.4$, $0.8$, $1.6$ and $3.2$ $pN$. Large stretching force $f$ results in more
compact plectonemes, i.e., reduces the radius of the plectonemic helix(Fig 4).
Indeed, in the coexistence state of DNA supercoiling, the equilibrium excess
linking number density of the P phase is
$\sigma_{P}^{eq}=\frac{1}{c_{P}}(\frac{-2c_{P}(f-\sqrt{\frac{k_{B}Tf}{l_{b}}})}{1-c_{P}/c_{S}})^{\frac{1}{2}}\sim\sqrt{f/l_{p}}$Marko
(2007), which suggests that $\sigma_{P}^{eq}$ increases with $f$.

At high salt concentration, the screening of the Coulomb repulsion potential
reduces the effective repulsive DNA diameter, leading to the shorter torsional
persistence length of the plectoneme $l_{p}$ and smaller radius of the
plectonemic helixRybenkov  _et al._ (1997). Meanwhile, based on the
coexistence state relation $\sigma_{P}^{eq}\propto\sqrt{f/l_{p}}$, lowering
$l_{p}$ increases the $\sigma_{P}$, i.e., making the plectonemes more compact
(Fig. 4). Accordingly, we adjusted the torsional persistence length of
plectoneme $l_{p}$ to represent the impacts from the ionic strengthRybenkov
_et al._ (1997), specifically, $l_{p}=25$, $23$ and $21$ $nm$ to represent the
low, regular and high salt concentrations, respectively.

According to Kramers’ reaction-rate theory, the nucleation rate is
$k_{nuc}\propto exp(-\beta\epsilon_{nuc})$ Daniels and Sethna (2011), where
the energy barrier $\epsilon_{nuc}\propto\sqrt{f}$, which incorporates the
bending energy of the loop and the loss of the elastic energy stored
originally in the stretched DNA(Fig. 4)Marko and Neukirch (2012). Therefore,
the large forces can lower the nucleation rate. Again, according to the
relation $\sigma_{P}^{eq}\propto\sqrt{f/l_{p}}$, large force and small $l_{p}$
(high salt concentration) have a similar effect on the nucleation rate.
Indeed, the nucleation rates obtained from simulations (Fig. 4) suggest that
large stretching forces $f$ and the high salt concentrations (small $l_{p}$)
can suppress the nucleation events. We derived the nucleation relation (Text
S8)

$k_{nuc}\propto exp(-2\sqrt{f}C(l_{p}^{\star}-l_{p})),$ (D1)

where $2C$ and $l_{p}^{\star}$ are fitting parameters (dashed curves in Fig.
4). As $l_{p}\rightarrow l_{p}^{\star}$ (salt concentration decreases), the
nucleation rate increases. That is because low salt concentration increases
the radius of the plectonemic helices i.e.,$R_{plect}\rightarrow R_{loop}$
(bending energy of DNA in plectonemes increases), then lower the energy
barrier from plectonemic helix to loop (reversed buckling transition, red
arrow in Fig. 4). By fitting, we obtained $l_{p}^{\star}=27nm$.

The results above indicate that the plectonemes frequently vanish and re-
emerge (via nucleation), which lead to the hopping. The surviving time of the
individual plectoneme from nucleation to vanishing is defined as plectoneme
life-time. Fig. 4 shows that the mean life-time $\tau_{plec}$ depends on
stretching force and salt concentration. The large $f$ or small $l_{p}$
increases the plectoneme life-time. As discussed above, large $f$ or small
$l_{p}$ lowers the plectonemes nucleation and also leads to more compact
plectonemes. The resulting compactness stabilizes existing plectonemes, which
means long mean life-time of the individual plectonemes. By obtaining the
balance between the plectoneme nucleation and vanishing, we can derive a mean
life time of the individual plectoneme (see Text S8)

$\tau_{plec}\propto(2l_{t}-l_{p})fexp(C\sqrt{f}(l_{p}^{\star}-l_{p})).$ (D2)

where $C$ and $l_{p}^{\star}$ are the same as those in Eq D1.

#### Plectoneme hopping and the associated energetic changes

An example of hopping $over10kbp$ within about $0.1s$ is shown in Fig.4.
Plectoneme hopping is a special case of the DNA conformational
rearrangementvan Loenhout  _et al._ (2012). As addressed above, plectoneme
diffusion happens comparatively slowly, only hundreds of bps per second
($\sim\sqrt{2Dt}$). Hopping , however, is a rapid process, potentially,
facilitating fast site-specific recombination and enhancer activated
expressionSheinin and Wang (2012). Since plectoneme hopping is related to
nucleation, the large forces inhibiting the nucleation thus repress the
hopping.

We also found that plectoneme hopping is associated with significant energetic
changes between the two phases. We identify the nucleation as the starting
point of a hopping event. Then we collected numerous simulations (5000) that
contain the hopping events and aligned the starting points of the hopping
events together for averaging the energetics. The energy change of the S phase
accumulates up to about $20k_{B}T$ before the nucleation, while that of the P
phase falls down to about $20k_{B}T$. These significant energetic changes
facilitate the abrupt nucleating-vanishing processes (Fig. 4). After the
formation of the new plectoneme by nucleation, the energetic changes rapidly
converge to $0$. The energetic offsets between the two phases thus makes the
total energy change smaller as thermal fluctuations, i.e.,$1k_{B}T$(inset of
Fig. 4). It should be noted that there is no persistent energetic injection
into the system. Thus the energy change must be compensated by the
complementary reversed processes (one of the multiple plectonemes vanishes,
corresponding to reversed buckling transition shown in Fig. 4).

## Methods

Figure 5: The variables and tentative loops in the two-phase dynamics. (a) The
phase boundaries $X_{\alpha}^{l/r}$ describe the phase-transformation between
the stretched phase and the $\alpha$-th plectoneme, which involve the
transverse displacement of the $\alpha$-th plectoneme described by
$z_{\alpha}$, the displacement of the $\alpha$-th section of the stretched
phase and the slithering motion of parallel segments of the plectoneme
described by $v_{slt,\alpha}$. (b) The algorithm and the integrations are
conducted in the discrete bins along the DNA axes. The tentative loops occupy
the bins.

### The DNA topology and in-extensibility constraints

The Brownian dynamics Eq.B2 is completed with addition of two constraints. The
first one is the topology of supercoiling and the second one is the in-
extensibility of the model DNA. The DNA topology constraint, i.e., $\Delta
Lk=const.$, gives the allocation,

$\Delta Lk=\Delta Lk_{S}+\Delta Lk_{P},$ (M1)

where $\Delta Lk_{S}$ and $\Delta Lk_{P}$ are the excess linking numbers of S
and P phases, respectively. They can be written as the functions of $X$ and
$\sigma_{S/P}^{eq}$, i.e., $\Delta
Lk_{S}=\sum_{\alpha=0}\frac{\Omega_{0}}{2\pi}\int_{X_{\alpha}^{r}}^{X_{\alpha+1}^{l}}\sigma_{S}^{eq}(u,\mathbf{X})du$
and $\Delta
Lk_{P}=\sum_{\alpha=0}\frac{\Omega_{0}}{2\pi}\int_{X_{\alpha+1}^{l}}^{X_{\alpha+1}^{r}}\sigma_{P}^{eq}(\mathbf{X})du$,
where $\sigma_{S}^{eq}$ and $\sigma_{P}^{eq}$ satisfy the torque equilibrium
$c_{S}\sigma_{S}^{eq}(u,\mathbf{X})=c_{P}\sigma_{P}^{eq}(X)=T_{q}^{in}$.

The second constraint originates from the in-extensibility of the model DNA.
Eq B2 describes the phase-transformation motion perpendicular to the
stretching direction. As shown in Fig. 5, $X_{\alpha}^{l,r}$ specify the
$\alpha$-th plectoneme, and also involve the displacements of $\alpha$-th and
$\alpha+1$-th sections of the stretched phase. This connectivity between the
plectonemic and stretched phases additionally confines the translational
displacements parallel to the stretching direction. Plectonemes and stretched
phase are subjected to the tensions at boundaries
$f_{u}|_{u=X_{\alpha}^{l/r}}$, on which potential $\Phi_{0}(\mathbf{X})$
depends. The components of motion parallel to the stretching direction,
involve the linear combinations of $X_{\alpha}^{l/r}$. Notice that the
coordinate of the $\alpha$-th plectoneme is
$z_{\alpha}=X_{\alpha}^{l}-\sum_{1}^{\alpha-1}l_{plect,\alpha’}$, where
$l_{plect,\alpha’}=X_{\alpha’}^{r}-X_{\alpha’}^{l}$ is the contour length of
the $\alpha’$-th plectoneme. The velocity of the $\alpha$-th stretched section
is defined as $\dot{z}_{\alpha}-\dot{X}_{\alpha}^{l}$. The internal slithering
velocity of the parallel segments of the plectoneme
$v_{slt,\alpha}=\frac{1}{2}(\dot{X}_{\alpha}^{l}+\dot{X}_{\alpha}^{r})$. The
constraints on parallel components result in mechanical balances on the
stretching direction

$\begin{split}&f_{u}|_{u=X_{\alpha}^{r}}-f_{u}|_{u=X_{\alpha}^{l}}-\mu_{slt}v_{slt,\alpha}=\gamma_{2,\alpha}\dot{z}_{\alpha}+\sqrt{2\beta^{-1}\gamma_{2,\alpha}}\dot{W}_{3}\\\
&f_{u}|_{u=X_{\alpha}^{l}}-f_{u}|_{u=X_{\alpha-1}^{r}}=\gamma_{3,\alpha}(\dot{z}_{\alpha}-\dot{X}_{\alpha}^{l})+\sqrt{2\beta^{-1}\gamma_{3,\alpha}}\dot{W}_{4}\end{split},$
(M2)

where $\mu_{slt}$ is slithering drag coefficient and $\gamma_{2,\alpha}$ is
the drag coefficient for the transverse displacement of the $\alpha$-th
plectoneme, $\dot{W}_{3}$ is the white noises on the displacement,
$\gamma_{3,\alpha}$ the drag coefficient for translational displacement of
$\alpha$-th stretched section, $\dot{W}_{4}$ is the white noises on the motion
of the $\alpha$-th stretched section. The first relation in Eq M2 describes
the constraints on the transverse displacement of the $\alpha$-th plectoneme,
and the second one describes constraints between the $\alpha-1$ and
$\alpha$-th plectonemes. In the dWLC method, the effective viscous drags is
$6\pi\eta R_{H}$ for each segment. Thus $\gamma_{1,\alpha}=6\pi\eta
R_{H}\frac{l_{plect,\alpha}}{l_{0}}\times\frac{1}{2}=0.6\pi\eta
l_{plect,\alpha}$ where we have used $\frac{R_{H}}{l_{0}}\simeq 0.2$ and
$\frac{1}{2}$ means a half of plectoneme contour length in Eq B2,
$\gamma_{2,\alpha}=1.2\pi\eta l_{plect,\alpha}$ and
$\gamma_{3,\alpha}=1.2\pi\eta(X_{\alpha}^{l}-X_{\alpha-1}^{r})$.

Neverthless, the friction coefficient $\mu_{slt}$ for slithering is unknown.
According to a relation of the slithering drag,
$\mu_{slt}\propto\frac{l_{plect,\alpha}}{ln(\frac{R_{plect}}{R})}$Marko and
Siggia (1994), where the plectoneme radius
$R_{plect}\approx\chi_{1}ln(\frac{\chi_{2}}{f})$ and $\chi_{1}=0.4nm$ and
$\chi_{2}=360pN$van Loenhout  _et al._ (2012). Calibration of the two-phase
model with the dWLC method gives
$\mu_{slt}=\frac{0.4l_{plect,\alpha}}{ln(\frac{R_{plect}}{R})}$. For
supercoiling with a single plectoneme, on average, $\dot{z}=v_{slt}$, and then
the viscous drag is $\gamma_{1,1}+\gamma_{2,1}+\mu_{slt,1}$. According to the
Einstein relation, the diffusion coefficient of individual plectoneme is
$D=\frac{k_{B}T}{\gamma_{1,1}+\gamma_{2,1}+\mu_{slt,1}}$.

### The routine for simulating the two-phase dynamics

$i$. Generating the conformation $\\{X_{\alpha}^{l/r},\sigma_{S/P}^{eq}\\}$.

The model DNA is discretized as bins with a size $5nm$, on which the
integrations are numerically conducted (Fig. 5). The observable $\mathbf{X}$
follows Eq B2 under constraints Eqs M1 and M2. We select the time-step as
$dt=2.5\frac{\zeta_{R}L^{2}}{k_{B}Tl_{t}}$, specifically, $\Delta
t=0.1,0.2,0.4,0.9$ and $5ms$ for $3,4.5,6,9$ and $21kb$, respectively. At each
time-step, Eq B2 and constraints Eqs M1 and M2 together constitute the self-
consistent equations that can be numerically solved with iteration procedure.
Take 20 iteration cycles at each time-step, and the convergence (relative
deviation of $\sigma_{S/P}<0.001$) is achieved within 10 iteration cycles.

$ii$. Sampling new plectonemes. The energy of the supercoiled DNA with n
tentative loops for nucleation is

$\Phi_{n}(\mathbf{X})=\Phi_{0}(\mathbf{X})+\sum_{i}^{n}E_{loop,i},$ (M3)

where $E_{loop,i}=8k_{B}T\sqrt{\frac{l_{b}f_{u}}{k_{B}T}}$ is the energy
penalty for looping at location $u$Marko and Neukirch (2012). The excess
linking number of the $n$ tentative loops is $\Delta
Lk_{l}=n+\frac{1}{2\pi}\frac{\Delta\theta}{\Delta u}\sum L_{loop,i}$, where
the relative twist is proportional to torque, i.e.,
$\frac{\Delta\theta}{\Delta u}=\frac{T_{q}^{in}}{k_{B}Tl_{t}}$. The time-scale
for looping is on microsecondsDaniels and Sethna (2011), much shorter than
$dt$. We thus assume that a portion of $\Delta Lk_{S}$ is stored in loops,
then the tentative state follows the two constraints $L_{S}=L_{S}’+\sum
L_{loop,i}$ and $\Delta Lk_{S}=\Delta Lk_{S}’+\Delta Lk_{l}$ where $L_{S}’$
and $Lk_{S}’$ are the contour length and linking number of the S phase of the
tentative state, respectively, and the contour length of a loop is
$L_{loop,i}=\sqrt{2\pi^{2}}\sqrt{k_{B}Tl_{b}/f_{u}}$ Marko and Neukirch
(2012).

Figure 6: Flow chart of simulation routine.

The loops along the stretched DNA can be regarded as hard-core particles
occupying a 1-dimensional discrete bins (Fig. 5), where each bin characterizes
the size of the thermal fluctuations of the DNA. The length of each bin
$d_{bin}\approx 5nm$, defined as the transverse amplitude of the
fluctuationsMarko and Siggia (1995b); Marko and Neukirch (2012). Define a
partition function of the tentative state with n loops,

$\mathcal{P}_{n}=\sum_{I_{1}<I_{2}<\cdot\cdot\cdot<I_{n}}exp[-\beta\Phi_{n}(\mathbf{X})],$
(M4)

where $I_{i}$ denotes the location of the i-th tentative loop. Accordingly,
$P_{0}=\frac{\mathcal{P}_{0}}{\sum_{i}\mathcal{P}_{i}}$ is the probability of
the supercoiling without tentative loops and
$P_{n}=\frac{\mathcal{P}_{n}}{\sum_{i}\mathcal{P}_{i}}$ is the probability of
supercoiling with $n$ tentative loops. According to the kinetic Monte Carlo
methodVan Santen and Sautet (2015), generate a random $\mathcal{R}$ in (0,1).
The range (0,1) is divided into $(0,P_{0}]$,
$(P_{0},P_{0}+P_{1}]$,$\cdot\cdot\cdot$,$(\sum_{i}^{\mathcal{N}-1}P_{i},\sum_{i}^{\mathcal{N}}P_{i}],\cdot\cdot\cdot$,
where $\sum_{\mathcal{N}+1}^{\infty}P_{i}\simeq 0$ and
$\sum_{i}^{\mathcal{N}}P_{i}\simeq 1$ (often $\mathcal{N}\leqslant 2$). If
$\mathcal{R}$ lies in $(\sum_{i}^{N-1}P_{i},\sum_{i}^{N}P_{i}]$, the state
with n loops is sampled. The same procedure is applied to locate the $n$ loops
according to the distribution
$\frac{exp[-\beta\Phi_{n}(\mathbf{X})]}{\mathcal{P}_{n}}$. (For a fully
stretched DNA, distribution Eq.M4 implies an entropy brought by $n$ loops,
i.e., $S_{n}=k_{B}nln(\frac{L-nL_{loop}}{d_{bin}})-k_{B}ln(n!)$, the
$\frac{L-nL_{loop}}{d_{bin}}$ bins occupied by $n$ identical particles Marko
and Neukirch (2012).)

$iii$. Removing plectonemes

If $\mathcal{N}’$ plectonemes satisfy
$l_{plect,\alpha}<\sqrt{2\pi^{2}}\sqrt{k_{B}Tf/l_{b}}$ (often
$\mathcal{N}’\leqslant 2$). Similar to $ii$, remove the plectonemes according
to the distribution $exp[-\beta\Phi_{n}(\mathbf{X})]$. $ii$ and $iii$ generate
and remove plectonemes according to the Boltzmann distribution.

Update the conformation $\\{X_{\alpha}^{l/r},\sigma_{S/P}^{eq}\\}$, and then
repeat $i$, $ii$ and $iii$. All the steps can be organized as the flow chart
(Fig. 6).

## Discussion and Conclusion

In this study, based on the polymer physics or the WLC model of DNA, we have
constructed a two-phase dynamic model of the DNA supercoiling. To establish
such a model, we demonstrate the fast dynamics of torque equilibrium within
the S and P (plectonemic) phases, identify the phase boundaries as slow
observables, and derive the energetics associated with the slow observables.
By comparing with the WLC model, we have shown that the two-phase dynamic
model provides a physically justified and efficient way in representing DNA
supercoiling dynamics. The model is particularly suited for studying the
plectoneme dynamics from milliseconds to seconds. Correspondingly, we have
demonstrated the plectoneme diffusion, nucleation, and hopping, consistently
with measurements from single molecule experimentsShashkova and Leake (2017);
Sheinin and Wang (2012); van Loenhout  _et al._ (2012).

In order to probe the time-scale separation, we have investigated the dynamics
of torque approaching to equilibrium within the pure S and P phases,
respectively. Based on the WLC model, the internal torque transports as twist
angle propagates. Similarly, by introducing an effective twist angle of the
plectonemic helix, the internal torque propagation within the P phase can be
described, comparably. Such derivations are confirmed by numerical simulations
based on the dWLC, which consistently show that the torque transport is fast,
i.e., on either the S or P phase of the supercoiled DNA.

Accordingly, we choose the the boundaries of the S and P phases as observables
to describe the slow dynamics. Based on the interphase torque equilibrium, one
can then derive the free energy associated with the collective phase
boundaries, and the corresponding equation of motion via the Langevin
dynamics.

To compare torsional responses from the two-phase model with that from the WLC
model, we tested at a constant rate of supercoil accumulation, $\Delta Lk=10$
per second, comparable to the RNA polymerization rate during transcription
elongation (i.e.,$\sim$ 100 bp per sec), in which $\Delta Lk=0.1$ is injected
into DNA for the RNA polymerase enzyme unwinding to synthesize one
nucleotideThomen  _et al._ (2005); Jing  _et al._ (2018). The discontinuities
on torque and extension curves are reproduced by the two-phase dynamic model,
which reflect the buckling transition during plectoneme nucleationForth  _et
al._ (2008); Marko and Neukirch (2012); Emanuel  _et al._ (2013); Daniels and
Sethna (2011); Walker  _et al._ (2018); Ott  _et al._ (2020). The torque and
extension versus linking number curves suggest that the non-equilibrium effect
of the supercoiling accumulation (from the growth of the plectonemic phase)
induced by RNAP is constantly relaxed during each polymerization step or
elongation cycle. Hence, such an RNA elongation process can be regarded as
quasi-steady state. An obvious advantage of the currently developed two-phase
dynamic model of DNA supercoiling comes from its trivial computational cost as
we only deal with the slow degrees of freedom of the system. Besides, we
performed simulations of the plectonemes growth under a constant torque since
some studies also regarded RNAPs as a torsional motor which imposes a constant
torqueMielke  _et al._ (2004). We accordingly determined frictional
coefficient due to slithering in plectoneme by calibrating with the WLC model.

With the two-phase dynamic model we described the DNA supercoiling dynamics
and particularly reproduced the plectoneme diffusion, nucleation and hopping.
The plectoneme diffusion coefficient is about 0.1 to 0.2$\mu m^{2}/s$,
decreasing slightly with the stretching force. Large force indeed makes the
plectoneme more compact, which increases the friction for the slithering of
parallel segments of the plecoenemes. Consequently, the large force reduces
the diffusive capability of plectoneme. Large stretching force also suppresses
the nucleation rate, as the force increases the energy penalty for forming a
loop. The ability of salt concentration to screen the electrostatic repulsion
between opposing segments in plectonemes also affects the radius of
plectonemic helixRybenkov  _et al._ (1997); van Loenhout  _et al._ (2012).
Similar to large stretching force, strong ionic strength results in compact
plectonemes. The resulting compact plectonemes require higher energy penalty
to fluctuate, i.e., large stretching force or strong ionic strength can thus
stabilize plectonemes. Accordingly, large force or strong ionic strength can
prolong the mean life-time of the individual plectoneme. We also observed that
a plectoneme can perform a long distance hopping. Specifically, the abruptly
conformational change is accompanied by significant energetic changes.
However, the total energetic change of the two phase is trivial or as small as
thermal fluctuations.

Recent experiments reveal the long-distance cooperative and antagonistic
dynamics between multiple RNAPsKim  _et al._ (2019). The experimental
evidences suggest that the superpcoil cancellation between adjacent RNAPs
makes similar apparent transcription elongation rate. However, upon promoter
repression or divergently transcribed genes, RNAPs switch their behavior from
collaborative to antagonistic due to the build-up supercoils. There are two
important dynamic behaviors that are incorporated in our two-phase model.
First, RNAPs can be sensitive to the torque close to the stalling condition
(5$\sim$11 $pN$ $nm$)Ma  _et al._ (2013, 2019). A static coexistence torque
may not slow down or even stall the transcription elongation, but the torque
fluctuation (time-dependent) may. Thus time-dependent torque may induce
effects on the long-distance cooperative behaviors of RNAPs. Second, supercoil
can fast propagate as twist (torque transport) or slowly propagate as writhe
(e.g, plectoneme diffusion). These two kinds of supercoiling propagations are
often confounded in modeling supercoiling-dependent dynamics. During
transcription elongation, the build-up supercoil can be stored in writhe and
can be transmitted as twist that can be felt by RNAPs as the torsional
stressMa  _et al._ (2019); Ma and Wang (2014b, a); Ma  _et al._ (2013). The
two-phase dynamic model of DNA supercoiling is expected to be applied to
improve understanding the supercoiling-dependent transcription.

DNA supercoiling facilitates the site juxtaposition of two distal DNA sites by
slithering opposing segments of the interwound superhelix or bending DNA to
coilOram  _et al._ (1997); Polikanov  _et al._ (2007); Parker and Halford
(1991); Marko (1997b), it then promotes protein binding or enhance the gene
expressionvan den Broek  _et al._ (2008); Liu  _et al._ (2001); Lomholt  _et
al._ (2009). As we have seen, slithering motion of the segments in a
plectomeic helix is accompanied with plectoneme diffusion. However, the
diffusion may not offer an efficient mechanism for site juxtaposition. An
efficient mechanism for achieving site juxtaposition is combining plectoneme
diffusion and hopping. These mechanisms can be efficiently explored by
employing the two phase dynamics of the DNA supercoiling which is physical
justified and computationally trivial.

## Acknowledgements

This work has been supported by NSFC Grant $\\#11775016$ and $\\#11635002$. JY
has been supported by the CMCF of UCI via NSF DMS 1763272 and the Simons
Foundation grant $\\#594598$ and start-up fund from UCI. We acknowledge the
computational support from the Beijing Computational Science Research Center
(CSRC). Thanks Dr Xinliang Xu for discussions on the work.

## References

  * Bates  _et al._ (2005) A. D. Bates, A. Maxwell, T. Maxwell, _et al._ ,  _DNA topology_ (Oxford University Press, USA, 2005).
  * Fuller (1971) F. B. Fuller, Proceedings of the National Academy of Sciences 68, 815 (1971).
  * Koster  _et al._ (2010) D. A. Koster, A. Crut, S. Shuman, M.-A. Bjornsti, and N. H. Dekker, Cell 142, 519 (2010).
  * Lavelle (2014) C. Lavelle, Current opinion in genetics & development 25, 74 (2014).
  * Dixon  _et al._ (2012) J. R. Dixon, S. Selvaraj, F. Yue, A. Kim, Y. Li, Y. Shen, M. Hu, J. S. Liu, and B. Ren, Nature 485, 376 (2012).
  * Glinsky (2015) G. Glinsky, arXiv preprint arXiv:1507.05368 (2015).
  * Stadhouders  _et al._ (2019) R. Stadhouders, G. J. Filion, and T. Graf, Nature 569, 345 (2019).
  * Naughton  _et al._ (2013) C. Naughton, N. Avlonitis, S. Corless, J. G. Prendergast, I. K. Mati, P. P. Eijk, S. L. Cockroft, M. Bradley, B. Ylstra, and N. Gilbert, Nature structural & molecular biology 20, 387 (2013).
  * Gilbert and Allan (2014) N. Gilbert and J. Allan, Current opinion in genetics & development 25, 15 (2014).
  * Corless and Gilbert (2016) S. Corless and N. Gilbert, Biophysical reviews 8, 245 (2016).
  * Postow  _et al._ (2004) L. Postow, C. D. Hardy, J. Arsuaga, and N. R. Cozzarelli, Genes & development 18, 1766 (2004).
  * Liu and Wang (1987) L. F. Liu and J. C. Wang, Proceedings of the National Academy of Sciences 84, 7024 (1987).
  * Ma  _et al._ (2013) J. Ma, L. Bai, and M. D. Wang, Science 340, 1580 (2013).
  * Chong  _et al._ (2014) S. Chong, C. Chen, H. Ge, and X. S. Xie, Cell 158, 314 (2014).
  * Burnier  _et al._ (2008) Y. Burnier, J. Dorier, and A. Stasiak, Nucleic acids research 36, 4956 (2008).
  * Schvartzman  _et al._ (2013) J. B. Schvartzman, M.-L. Martínez-Robles, P. Hernández, and D. B. Krimer, Biochemical Society Transactions 41, 646 (2013).
  * Hirano (2000) T. Hirano, Annual review of biochemistry 69, 115 (2000).
  * Strick  _et al._ (2004) T. R. Strick, T. Kawaguchi, and T. Hirano, Current biology 14, 874 (2004).
  * Kim  _et al._ (2019) S. Kim, B. Beltran, I. Irnov, and C. Jacobs-Wagner, Cell 179, 106 (2019).
  * Teves and Henikoff (2014) S. S. Teves and S. Henikoff, Nucleus 5, 211 (2014).
  * Tantale  _et al._ (2016) K. Tantale, F. Mueller, A. Kozulic-Pirher, A. Lesne, J.-M. Victor, M.-C. Robert, S. Capozi, R. Chouaib, V. Bäcker, J. Mateos-Langerak, _et al._ , Nature communications 7, 1 (2016).
  * Jing  _et al._ (2018) X. Jing, P. Loskot, and J. Yu, Physical biology 15, 066007 (2018).
  * Kimura and Hirano (1997) K. Kimura and T. Hirano, Cell 90, 625 (1997).
  * Swedlow and Hirano (2003) J. R. Swedlow and T. Hirano, Molecular cell 11, 557 (2003).
  * Seol and Neuman (2016) Y. Seol and K. C. Neuman, Biophysical reviews 8, 101 (2016).
  * Ma and Wang (2014a) J. Ma and M. Wang, Transcription 5, e28636 (2014a).
  * Ma and Wang (2014b) J. Ma and M. D. Wang, Cell Cycle 13, 337 (2014b).
  * Fuller (1978) F. B. Fuller, Proceedings of the National Academy of Sciences 75, 3557 (1978).
  * Parker and Halford (1991) C. N. Parker and S. E. Halford, Cell 66, 781 (1991).
  * Liu  _et al._ (2001) Y. Liu, V. Bondarenko, A. Ninfa, and V. M. Studitsky, Proceedings of the National Academy of Sciences 98, 14883 (2001).
  * Huang  _et al._ (2001) J. Huang, T. Schlick, and A. Vologodskii, Proceedings of the National Academy of Sciences 98, 968 (2001).
  * Polikanov  _et al._ (2007) Y. S. Polikanov, V. A. Bondarenko, V. Tchernaenko, Y. I. Jiang, L. C. Lutter, A. Vologodskii, and V. M. Studitsky, Biophysical journal 93, 2726 (2007).
  * Marko (1997a) J. F. Marko, Physical Review E 55, 1758 (1997a).
  * Strick  _et al._ (1999) T. Strick, J.-F. Allemand, D. Bensimon, R. Lavery, and V. Croquette, Physica A: Statistical Mechanics and its Applications 263, 392 (1999).
  * Ghatak and Mahadevan (2005) A. Ghatak and L. Mahadevan, Physical review letters 95, 057801 (2005).
  * Marko (2007) J. F. Marko, Physical Review E 76, 021926 (2007).
  * Forte  _et al._ (2019) G. Forte, M. Caraglio, D. Marenduzzo, and E. Orlandini, Physical Review E 99, 052503 (2019).
  * Marko and Siggia (1994) J. F. Marko and E. D. Siggia, Science 265, 506 (1994).
  * Bustamante  _et al._ (1994) C. Bustamante, J. Marko, E. Siggia, and S. Smith, Science 265, 1599 (1994).
  * Strick  _et al._ (1996) T. R. Strick, J.-F. Allemand, D. Bensimon, A. Bensimon, and V. Croquette, Science 271, 1835 (1996).
  * Strick  _et al._ (1998) T. Strick, J.-F. Allemand, D. Bensimon, and V. Croquette, Biophysical journal 74, 2016 (1998).
  * Bustamante  _et al._ (2000) C. Bustamante, S. B. Smith, J. Liphardt, and D. Smith, Current opinion in structural biology 10, 279 (2000).
  * Strick  _et al._ (2000) T. Strick, J.-F. Allemand, V. Croquette, and D. Bensimon, Progress in biophysics and molecular biology 74, 115 (2000).
  * Koster  _et al._ (2005) D. A. Koster, V. Croquette, C. Dekker, S. Shuman, and N. H. Dekker, Nature 434, 671 (2005).
  * van Loenhout  _et al._ (2012) M. T. van Loenhout, M. De Grunt, and C. Dekker, Science 338, 94 (2012).
  * Bustamante  _et al._ (2003) C. Bustamante, Z. Bryant, and S. B. Smith, Nature 421, 423 (2003).
  * Forth  _et al._ (2008) S. Forth, C. Deufel, M. Y. Sheinin, B. Daniels, J. P. Sethna, and M. D. Wang, Physical review letters 100, 148301 (2008).
  * Marko and Siggia (1995a) J. F. Marko and E. D. Siggia, Macromolecules 28, 8759 (1995a).
  * Marko and Siggia (1995b) J. Marko and E. Siggia, Physical Review E 52, 2912 (1995b).
  * Moroz and Nelson (1997) J. D. Moroz and P. Nelson, Proceedings of the National Academy of Sciences 94, 14418 (1997).
  * Phillips  _et al._ (2009) R. R. B. Phillips, J. Kondev, J. Theriot, _et al._ ,  _Physical biology of the cell_ (Garland Science, 2009).
  * Marko (1998) J. F. Marko, Physical Review E 57, 2134 (1998).
  * Daniels and Sethna (2011) B. C. Daniels and J. P. Sethna, Physical Review E 83, 041924 (2011).
  * Marko and Neukirch (2012) J. F. Marko and S. Neukirch, Physical Review E 85, 011908 (2012).
  * Emanuel  _et al._ (2013) M. Emanuel, G. Lanzani, and H. Schiessel, Physical Review E 88, 022706 (2013).
  * Mitchell and Harris (2013) J. S. Mitchell and S. A. Harris, Physical review letters 110, 148105 (2013).
  * Mitchell  _et al._ (2011) J. Mitchell, C. Laughton, and S. A. Harris, Nucleic acids research 39, 3928 (2011).
  * Irobalieva  _et al._ (2015) R. N. Irobalieva, J. M. Fogg, D. J. Catanese, T. Sutthibutpong, M. Chen, A. K. Barker, S. J. Ludtke, S. A. Harris, M. F. Schmid, W. Chiu, _et al._ , Nature communications 6, 1 (2015).
  * Matek  _et al._ (2015) C. Matek, T. E. Ouldridge, J. P. Doye, and A. A. Louis, Scientific reports 5, 7655 (2015).
  * Sutthibutpong  _et al._ (2016) T. Sutthibutpong, C. Matek, C. Benham, G. G. Slade, A. Noy, C. Laughton, J. P. K. Doye, A. A. Louis, and S. A. Harris, Nucleic acids research 44, 9121 (2016).
  * Hirsh  _et al._ (2013a) A. D. Hirsh, M. Taranova, T. A. Lionberger, T. D. Lillian, I. Andricioaei, and N. Perkins, Biophysical journal 104, 2058 (2013a).
  * Dans  _et al._ (2016) P. D. Dans, J. Walther, H. Gómez, and M. Orozco, Current opinion in structural biology 37, 29 (2016).
  * Goyal  _et al._ (2005) S. Goyal, N. C. Perkins, and C. L. Lee, Journal of Computational Physics 209, 371 (2005).
  * Lillian  _et al._ (2011) T. D. Lillian, M. Taranova, J. Wereszczynski, I. Andricioaei, and N. Perkins, Biophysical journal 100, 2016 (2011).
  * Hirsh  _et al._ (2013b) A. D. Hirsh, T. D. Lillian, T. A. Lionberger, M. Taranova, I. Andricioaei, and N. Perkins, Journal of Computational and Nonlinear Dynamics 8 (2013b).
  * Vologodskii (2006) A. Vologodskii, in  _Computational Studies of RNA and DNA_ (Springer, 2006) pp. 579–604.
  * Vologodskii  _et al._ (1992) A. V. Vologodskii, S. D. Levene, K. V. Klenin, M. Frank-Kamenetskii, and N. R. Cozzarelli, Journal of molecular biology 227, 1224 (1992).
  * Klenin  _et al._ (1991) K. V. Klenin, A. V. Vologodskii, V. V. Anshelevich, A. M. Dykhne, and M. D. Frank-Kamenetskii, Journal of molecular biology 217, 413 (1991).
  * Klenin  _et al._ (1998) K. Klenin, H. Merlitz, and J. Langowski, Biophysical journal 74, 780 (1998).
  * Bell  _et al._ (2012) D. Bell, P. Selokar, and T. D. Lillian, Biophysical Journal 102, 276a (2012).
  * Babamohammadi and Lillian (2020) S. Babamohammadi and T. D. Lillian, Biophysical Journal 119, 2517 (2020).
  * Walker  _et al._ (2018) P. U. Walker, W. Vanderlinden, and J. Lipfert, Physical Review E 98, 042412 (2018).
  * Ott  _et al._ (2020) K. Ott, L. Martini, J. Lipfert, and U. Gerland, Biophysical journal 118, 1690 (2020).
  * Wocjan  _et al._ (2009) T. Wocjan, K. Klenin, and J. Langowski, The journal of physical chemistry B 113, 2639 (2009).
  * Ivenso and Lillian (2016) I. D. Ivenso and T. D. Lillian, Biophysical journal 110, 2176 (2016).
  * Wada and Netz (2009) H. Wada and R. R. Netz, EPL (Europhysics Letters) 87, 38001 (2009).
  * Fasnacht  _et al._ (2004) M. Fasnacht, R. H. Swendsen, and J. M. Rosenberg, Physical Review E 69, 056704 (2004).
  * Zhang and Ma (2010) C. Zhang and J. Ma, The Journal of chemical physics 132, 244101 (2010).
  * Crut  _et al._ (2007) A. Crut, D. A. Koster, R. Seidel, C. H. Wiggins, and N. H. Dekker, Proceedings of the National Academy of Sciences 104, 11957 (2007).
  * Rybenkov  _et al._ (1997) V. V. Rybenkov, A. V. Vologodskii, and N. R. Cozzarelli, Nucleic acids research 25, 1412 (1997).
  * Neukirch and Marko (2011) S. Neukirch and J. F. Marko, Physical review letters 106, 138104 (2011).
  * Weinan (2011) E. Weinan,  _Principles of multiscale modeling_ (Cambridge University Press, 2011).
  * Mattis and Swendsen (2008) D. C. Mattis and R. H. Swendsen,  _Statistical Mechanics Made Simple_ (World Scientific Publishing Company, 2008).
  * Thomen  _et al._ (2005) P. Thomen, P. J. Lopez, and F. Heslot, Physical Review Letters 94, 128102 (2005).
  * Mielke  _et al._ (2004) S. P. Mielke, W. H. Fink, V. V. Krishnan, N. Grønbech-Jensen, and C. J. Benham, The Journal of chemical physics 121, 8104 (2004).
  * Sheinin and Wang (2012) M. Y. Sheinin and M. D. Wang, Science 338, 56 (2012).
  * Van Santen and Sautet (2015) R. A. Van Santen and P. Sautet,  _Computational methods in catalysis and materials science: an introduction for scientists and engineers_ (John Wiley & Sons, 2015).
  * Shashkova and Leake (2017) S. Shashkova and M. C. Leake, Bioscience reports 37, BSR20170031 (2017).
  * Ma  _et al._ (2019) J. Ma, C. Tan, X. Gao, R. M. Fulbright, J. W. Roberts, and M. D. Wang, Proceedings of the National Academy of Sciences 116, 2583 (2019).
  * Oram  _et al._ (1997) M. Oram, J. F. Marko, and S. E. Halford, Journal of molecular biology 270, 396 (1997).
  * Marko (1997b) J. F. Marko, Physica A: Statistical Mechanics and its Applications 244, 263 (1997b).
  * van den Broek  _et al._ (2008) B. van den Broek, M. A. Lomholt, S.-M. Kalisch, R. Metzler, and G. J. Wuite, Proceedings of the National Academy of Sciences 105, 15738 (2008).
  * Lomholt  _et al._ (2009) M. A. Lomholt, B. van den Broek, S.-M. J. Kalisch, G. J. Wuite, and R. Metzler, Proceedings of the National Academy of Sciences 106, 8204 (2009).