Twist-bend coupling, twist waves and the shape of DNA loops

TL;DR

This paper combines analytical and numerical methods to study the shape and dynamics of short DNA loops, revealing twist oscillations caused by twist-bend coupling due to DNA's groove asymmetry, and demonstrates their robustness against thermal fluctuations.

Contribution

It introduces a simple analytical variational shape that accurately models DNA loops and elucidates the role of twist-bend coupling in their shape and dynamics.

Findings

Twist oscillations are observed in DNA loops due to twist-bend coupling.

The analytical model reproduces loop energy with high accuracy.

Twist waves are stable against thermal fluctuations and exhibit diffusive motion.

Abstract

By combining analytical and numerical calculations, we investigate the minimal-energy shape of short DNA loops of approximately $100$ base pairs (bp). We show that in these loops the excess twist density oscillates as a response to an imposed bending stress, as recently found in DNA minicircles and observed in nucleosomal DNA. These twist oscillations, here referred to as twist waves, are due to the coupling between twist and bending deformations, which in turn originates from the asymmetry between DNA major and minor grooves. We introduce a simple analytical variational shape, that reproduces the exact loop energy up to the fourth significant digit, and is in very good agreement with shapes obtained from coarse-grained simulations. We, finally, analyze the loop dynamics at room temperature, and show that the twist waves are robust against thermal fluctuations. They perform a normal diffusive motion, whose origin is briefly discussed.

Figures

Tables

Table 1. Table 1: Comparison between harmonic-loop energies and the exact values for three different “loop” geometries with varying distance d 𝑑 d between the endpoints. d = 0 𝑑 0 d=0 corresponds to the closed loops of Fig.  1 .

d		$c_1$	$c_3$	$\beta EL/l_{\text{b}}$	Error(%)
$0$	1st	2.4048	–	14.2694	1.53
	3rd	2.3703	-0.2808	14.0572	0.02
	Exact			14.0550
$L/10$	1st	2.2187	–	12.1459	0.98
	3rd	2.1977	-0.2135	12.0295	0.02
	Exact			12.0286
$L/5$	1st	2.0415	–	10.2837	0.63
	3rd	2.0288	-0.1612	10.2198	0.01
	Exact	–	–	10.2194

Equations

$$\beta E=\frac{l_{\text{b}}}{2}{\int }_0^L\text{d}s{\left(\frac{\text{d}{\hat{\mathbf{e}}}_3}{\text{d}s}\right)}^2,$$

$$\frac{\text{d}{\hat{\mathbf{e}}}_i}{\text{d}s}=\left(\mathbf{\Omega }+{\omega }_0{\hat{\mathbf{e}}}_3\right)\times {\hat{\mathbf{e}}}_i,$$

$$\beta E=\frac{1}{2}{\int }_0^L\text{d}s\left(A_1{\Omega }_1^2+A_2{\Omega }_2^2+C{\Omega }_3^2+2G{\Omega }_2{\Omega }_3\right).$$

$${\kappa }^2\equiv {\left(\frac{\text{d}{\hat{\mathbf{e}}}_3}{\text{d}s}\right)}^2={\Omega }_1^2+{\Omega }_2^2,$$

$${\Omega }_1$$

$${\Omega }_2$$

$${\Omega }_3$$

$${\tilde{A}}_2=A_2\left(1-\frac{G^2}{A_{2}C}\right),$$

$$\frac{1}{l_{\text{b}}}=\frac{1}{2}\left(\frac{1}{A_1}+\frac{1}{{\tilde{A}}_2}\right),$$

$${\hat{\mathbf{e}}}_3=\cos \theta \hat{\mathbf{x}}+\sin \theta \hat{\mathbf{y}},$$

$$\frac{\text{d}\theta (0)}{\text{d}s}=\frac{\text{d}\theta (L)}{\text{d}s}=0.$$

$$\frac{\text{d}{\theta }^{(1)}}{\text{d}s}=\frac{\pi c_1}{L}\sin \left(\frac{\pi s}{L}\right),$$

$$\beta E_{\text{HL}}^{(1)}={\left(\frac{\pi c_1}{2}\right)}^2\frac{l_{\text{b}}}{L}=14.2694\frac{l_{\text{b}}}{L},$$

$$\frac{\text{d}{\theta }^{(3)}}{\text{d}s}=\frac{\pi }{L}\left[c_1\sin \left(\frac{\pi s}{L}\right)+c_3\sin \left(\frac{3\pi s}{L}\right)\right],$$

$$\beta E_{\text{HL}}^{(3)}=\frac{{\pi }^2\left(c_1^2+c_3^2\right)}{4}\frac{l_{\text{b}}}{L}=14.0572\frac{l_{\text{b}}}{L},$$

$${\kappa }_{\max }=\frac{\pi }{L}\left(c_1-c_3\right).$$

$${\kappa }_{\text{turn}}=\frac{\pi }{L}\left(c_1\cos \alpha -c_3\cos 3\alpha \right),$$

$$\frac{\Delta \kappa }{{\kappa }_{\text{max}}}\equiv \frac{{\kappa }_{\text{max}}-{\kappa }_{\text{turn}}}{{\kappa }_{\text{max}}}=1-\frac{c_1\cos \alpha -c_3\cos 3\alpha }{c_1-c_3}.$$

$${\Omega }_1^{(1)}$$

$${\Omega }_2^{(1)}$$

$${\Omega }_3^{(1)}$$

$${\Omega }_1^{(3)}$$

$${\Omega }_2^{(3)}$$

$${\Omega }_3^{(3)}$$

$$\beta E_{\text{HL}}^{(1)}={\left(\frac{\pi c_1}{2}\right)}^2\frac{l_{\text{b}}}{L}+\beta \Delta E(\varphi ).$$

$${\tilde{\kappa }}_{\max }\equiv {\sqrt{{\Omega }_1^2+{\Omega }_2^2}|}_{s=L/2},$$

$$\frac{2{\tilde{A}}_2}{A_1+{\tilde{A}}_2}\le \frac{{\tilde{\kappa }}_{\max }}{{\kappa }_{\max }}\le \frac{2A_1}{A_1+{\tilde{A}}_2},$$

$$\mathbf{R}\cdot \hat{\mathbf{x}}={\int }_0^{L/2}{\hat{\mathbf{e}}}_3\cdot \hat{\mathbf{x}}\text{d}s={\int }_0^{L/2}\cos [\theta (s)]\text{d}s=0.$$

$$\beta E={\int }_0^L\left(\frac{A{\dot{\theta }}^2}{2}-\mu \cos \theta \right)\text{d}s,$$

$$A\ddot{\theta }=\mu \sin \theta ,$$

$$\frac{A}{2}{\dot{\theta }}^2+\mu \cos \theta =\Gamma .$$

$$F\left(\frac{\pi -\theta (s)}{2},k\right)=F\left(\frac{\pi -\alpha }{2},k\right)-\sqrt{\frac{\mu }{A}}\frac{s}{k},$$

$$F(\varphi ,k)={\int }_0^{\varphi }\frac{d\omega }{\sqrt{1-k^2{\sin }^2\omega }}.$$

Full text

I Introduction

DNA often forms loops due to the action of proteins, which bind at two distant sites along its sequence and bring them in close contact with each other Alberts et al. (2002). DNA loops play an important role in many biological processes such as transcription, recombination and duplication. The loop length can range from 100 base pairs (bp), in the lac and gal operons of E. Coli Cournac and Plumbridge (2013), to tenths of thousands bp in more complex organisms Krivega and Dean (2012). For lengths much longer than the DNA persistence length, $l_{\text{b}}\approx 50$ nm ($150$ bp), entropic contributions dominate, leading to strong fluctuations in the shape, while in the opposite limit, namely on loops of length comparable to or smaller than $l_{\text{b}}$, the loop assumes approximately its minimal-energy shape. Several studies have focused on the latter regime, not only due to its biological relevance, but also because it allows to investigate the mechanics of highly-deformed DNA Balaeff et al. (1999); Kulić and Schiessel (2003); Sankararaman and Marko (2005); Zhang et al. (2006); Lee et al. (2010); Wilson et al. (2010); Cherstvy (2011); Vafabakhsh and Ha (2012); Le and Kim (2014); Chen et al. (2014); Mulligan et al. (2015).

Mechanical models used for short loops typically neglect sequence-dependent effects, and describe DNA as an isotropic continuous elastic rod. This description fails to account for a coupling interaction present in real DNA, that connects the bending and twisting degrees of freedom, and originates from the asymmetry of the DNA grooves Marko and Siggia (1994). The effect of such a twist-bend coupling interaction on the behavior of DNA at various scales has been recently discussed Nomidis et al. (2017); Skoruppa et al. (2017, 2018); Nomidis et al. (2019); Caraglio et al. (2019). The aim of this paper is to investigate its influence on the structure of short DNA loops. Starting from a simple wormlike chain (WLC) model, in which only the bending degrees of freedom are taken into account, we construct a simple variational ansatz for the loop shape, which we refer to as harmonic loop. This ansatz reproduces very accurately the exact loop shape, which is expressed in terms of elliptic integrals Yamakawa and Stockmayer (1972), while the exact loop energy is reproduced by the harmonic loop ansatz up to four significant digits. The advantage of the variational solution is that it involves simple trigonometric functions, from which various properties of the loop can be easily obtained. Combining with the results of previous work Skoruppa et al. (2018), we have extended the variational solution to more complex DNA models, with anisotropic bending and twist-bend coupling. The comparison with numerical simulations of various coarse-grained model of DNA shows that the harmonic-loop approximation performs extremely well in all cases. As such, it provides a basis for further analysis of the equilibrium and kinetic properties of the DNA loops, which are briefly addressed at the end of this paper.

II DNA elasticity and twist-bend coupling

The simplest continuum model of DNA is the WLC, in which a configuration is described by the tangent vector $\mathbf{\widehat{e}}_{3}(s)$, where $s$ is the curvilinear contour coordinate ($0\leq s\leq L$) and $L$ the total length. The energy takes the following form:

[TABLE]

where $\beta=1/k_{\text{B}}T$ is the inverse temperature and $l_{\text{b}}$ the bending persistence length. To include twist, the model can be extended by introducing two additional unit vectors $\mathbf{\widehat{e}}_{1}$ and $\mathbf{\widehat{e}}_{2}$, such that $\{\mathbf{\widehat{e}}_{1},\mathbf{\widehat{e}}_{2},\mathbf{\widehat{e}}_{3}\}$ forms an orthonormal basis. The vectors $\mathbf{\widehat{e}}_{1}$ and $\mathbf{\widehat{e}}_{2}$ keep track of the relative rotation around the $\mathbf{\widehat{e}}_{3}$ axis of neighboring points. In DNA, by convention, $\mathbf{\widehat{e}}_{2}$ connects the backbones of the two strands and $\mathbf{\widehat{e}}_{1}$ points towards the major groove. A generic configuration of the twisted rod can be then parametrized in terms of infinitesimal rotations connecting the orthonormal frame $\{\mathbf{\widehat{e}}_{1},\mathbf{\widehat{e}}_{2},\mathbf{\widehat{e}}_{3}\}$ in $s$ to a neighboring frame in $s+\text{d}s$. This rotation can be mathematically cast into the following differential equation

[TABLE]

where $i=1,2,3$, and $\omega_{0}\approx 1.75$ rad/nm is the intrinsic helical twist density. The solution of the previous equation for $\mathbf{\Omega}=\mathbf{0}$ corresponds to a twisted straight rod ($\mathbf{\widehat{e}}_{3}=\text{const}.$), in which $\mathbf{\widehat{e}}_{1}$ and $\mathbf{\widehat{e}}_{2}$ rotate with angular frequency $\omega_{0}$. One defines the three components of the vector $\mathbf{\Omega}$ along the frame as $\Omega_{i}\equiv\mathbf{\Omega}\cdot\mathbf{\widehat{e}}_{i}$. Here, $\Omega_{1}$ and $\Omega_{2}$ denote the bending densities along the two main axes of the molecule and $\Omega_{3}$ the excess twist density. From the analysis of the symmetry of a DNA molecule, Marko and Siggia derived the following continuum model Marko and Siggia (1994)

[TABLE]

The parameters $A_{1}$ and $A_{2}$ are the stiffnesses associated with bending over the backbone and the grooves, respectively, while $C$ is the intrinsic twist stiffness. Finally, the twist-bend coupling term $G$ leads to a correlation of the strain fields $\Omega_{2}$ and $\Omega_{3}$. Note that from Eq. (2) one obtains

[TABLE]

where $\kappa$ is the curvature. Using this relation, and setting $A_{1}=A_{2}=l_{\text{b}}$, while neglecting twist degrees of freedoms, one directly sees that model (3) reduces to (1).

The effect of twist-bend coupling on the mechanical properties of DNA has been investigated in a few recent papers Nomidis et al. (2017); Skoruppa et al. (2017, 2018); Caraglio et al. (2019); Nomidis et al. (2019). An interesting consequence of $G\neq 0$ is the existence of twist oscillations in curved DNA Skoruppa et al. (2018). In Ref. Skoruppa et al. (2018) the following minimal-energy shape of a minicircle was derived

[TABLE]

where $R$ is the average circle radius and

[TABLE]

an effective bending stiffness. Finally, $l_{\text{b}}$ is the bending persistence length, which within model (3) is given by

[TABLE]

i.e., the harmonic mean of the two stiffnesses Nomidis et al. (2017). The oscillations in the bending densities $\Omega_{1}$ and $\Omega_{2}$ arise from the geometrical constraints of the system, which in turn induce oscillations in the twist density $\Omega_{3}$ from the presence of twist-bend coupling ($G\neq 0$). These twist waves have been indeed observed in X-ray crystallographic structures of DNA bound to histone proteins Skoruppa et al. (2018). While Eq. (5) describes a torsionally-relaxed minicircle, this solution has also been recently extended to minicircles which are either over- or undertwisted Caraglio et al. (2019).

III Harmonic loops in the WLC model

Considerable attention has been devoted to the study of structural and dynamical properties of DNA loops in the past two decades Balaeff et al. (1999); Kulić and Schiessel (2003); Sankararaman and Marko (2005); Zhang et al. (2006); Lee et al. (2010); Wilson et al. (2010); Cherstvy (2011); Vafabakhsh and Ha (2012); Le and Kim (2014); Chen et al. (2014); Mulligan et al. (2015). Already in the early 70s, Yamakawa and Stockmayer Yamakawa and Stockmayer (1972) discussed the minimal-energy configuration of a semiflexible loop within the framework of the isotropic WLC [Eq. (1)]. The exact shape thus obtained is expressed in terms of inverse elliptic integrals, and its derivation is outlined in Appendix A. Although exact, these expressions are not easy to handle, therefore simpler approximate loop configurations have also been considered. For instance, in the context of DNA looping in the nucleosome, Kulic and Schiessel Kulić and Schiessel (2003) introduced a so-called circle-line approximation, in which the DNA conformation is built from straight segments and arcs of circles. A similar approach was followed by Sankararaman and Marko Sankararaman and Marko (2005), who studied DNA loops under tension. In the same spirit, we introduce here a different approximate shape for the loop, which we refer to as harmonic loop. Though still simple, it is found to be more accurate, and will allow us to construct a full analytical shape for the anisotropic model (3). From that, we can directly estimate several quantities of interest, such as twist oscillations, curvature variation and minimal energy.

As the problem is two-dimensional, one can describe the shape of the loop using a single parameter $\theta(s)$, defined as the angle the tangent $\mathbf{\widehat{e}}_{3}$ forms with the x-axis:

[TABLE]

where the unit vectors $\mathbf{\widehat{x}}$ and $\mathbf{\widehat{y}}$ lie on the plane of the loop. Since no rotational constraints are applied at the loop endpoints, the curvature $\kappa=|\text{d}\theta/\text{d}s|$ [see Eq. (4)] must vanish at the boundaries, hence

[TABLE]

A simple ansatz fulfilling these boundary conditions is

[TABLE]

which we refer to as first-order harmonic loop. Here, the dimensionless constant $c_{1}$ can be fixed by requiring that the endpoints coincide. This constraint is discussed in Appendix B and can be cast in the form $J_{0}(c_{1})=0$, with $J_{0}$ being the zeroth-order Bessel function of the first kind [see Eq. (30)]. The parameter $c_{1}$ is, thus, the first root of $J_{0}$, which is known to a high degree of accuracy, $c_{1}=2.40482556$.

Figure 1(a) shows with a red dashed line the loop shape, obtained by integrating Eq. (10), plugging into Eq. (8) and integrating once more. The red dashed line in Fig. 1(b) shows a plot of Eq. (10), multiplied by the loop length $L$, so as to render it dimensionless. For comparison, the same graphs show as blue solid lines the shapes and curvature corresponding to the exact solution (26) Yamakawa and Stockmayer (1972). There is a reasonable overall agreement between the exact solution and first-harmonic approximation. Plugging Eqs. (10) and (8) into Eq. (1), we find the following energy for the harmonic loop

[TABLE]

which is very close to the exact value $\beta E_{\text{exact}}=14.0550\,(l_{\text{b}}/L)$ of the Yamakawa-Stockmayer solution Yamakawa and Stockmayer (1972). The value in Eq. (11) improves upon the “circle-line” approximation of Ref. Kulić and Schiessel (2003), which gives $\beta E_{\text{CL}}=15.70\,(l_{\text{b}}/L)$.

Our approximation scheme can be systematically improved by including higher-order terms. We consider the following variational ansatz

[TABLE]

which extends Eq. (10) by adding an additional harmonic, consistent with the symmetry of the solution (we require $\text{d}\theta/\text{d}s$ to be symmetric around $s=L/2$, which excludes all even harmonics with frequency $2\pi n/L$ and $n$ integer). The two parameters $c_{1}$ and $c_{3}$ are now fixed by requiring both the closure of the loop and the energy minimization of Eq. (1). We, thus, find the values $c_{1}=2.3703$ and $c_{3}=-0.2808$ (more details can be found in Appendix B). The resulting shape and rescaled curvatures are shown as dashed red lines in Fig. 1(c) and (d), revealing an excellent agreement with the exact solution (solid blue line). The energy is found to be

[TABLE]

which matches the exact solution up to four significant digits. The ansatzes of Eqs. (10) and (12) can be extended to the case where the end-points are fixed at some finite distance $d$ (see Appendix B). Table 1 summarizes the optimal values of the coefficients $c_{1}$ and $c_{3}$ for the first- and third-harmonic approximations for some selected values of $d$. A comparison between the exact results and the harmonic approximations shows that the latter become even better with increasing $d$, i.e., as the distance between the end-points increases (for more details, see Appendix B).

From the simple analytical, but accurate, form of the loop shape one can easily obtain several estimates of the loop properties. Let us consider, for instance, the curvature $\kappa(s)=\text{d}\theta/\text{d}s$, which reaches its maximum value $\kappa_{\text{max}}$ at the apex $s=L/2$ of the loop. From the third-harmonic approximation (12), one finds

[TABLE]

Next, we wish to estimate the curvature $\kappa_{\text{turn}}$ at a distance of one helical repeat from the apex, i.e., at $s=L/2+2\pi/\omega_{0}$. Again, from Eq. (12) this is found to be

[TABLE]

where we introduced the angle $\alpha\equiv 2\pi^{2}/\omega_{0}L$. Finally, combining Eqs. (14) and (15), one finds the maximal curvature drop

[TABLE]

Consulting Table 1 for $d=0$ and considering a DNA loop of 100 bp ($L=34$ nm), yields $\kappa_{\text{max}}=0.24\,\text{nm}^{-1}$ and $\Delta\kappa/\kappa_{\text{max}}\approx 0.097$.

IV Minimal-energy configuration of DNA loops: Modulated twist waves

IV.1 Analytical results

The analysis of Eq. (16) for a loop of $100$ bp reveals a rather modest curvature drop at the scale of the helical-repeat length. This allows us to use Eq. (5), derived for a minicircle of average radius $R$, by replacing $1/R$ with the modulated harmonic-loop curvature. For instance, combining Eq. (10) with (5) yields the first-harmonic solution

[TABLE]

where a phase constant $\phi$ has been added, accounting for the torsional freedom of DNA at the boundaries (torsionally-unconstrained ends). Similarly, one can combine Eq. (5) with Eq. (12), so as to construct a more accurate approximation

[TABLE]

Similar to the minicircle case [Eqs. (5)], one notices the emergence of twist waves, originating from twist-bend coupling ($G\neq 0$). In this case, however, these are modulated by the varying curvature, which vanishes at the loop edges and is maximal at the loop apex.

By plugging Eqs. (17) into Eq. (3) and performing the integration in $s$, one can compute the total energy of the loop (see Appendix C for details)

[TABLE]

This expression is identical to Eq. (11) with the addition of a boundary term $\Delta E(\phi)$ depending on the phase $\phi$. One can show that this term is negligible for loops of about 100 bp, such as those considered here ($|\Delta E|/E^{(1)}_{\text{HL}}\sim 10^{-4}$). We, thus, conclude that the energy is quasi-degenerate, corresponding to an invariance of the double helix with respect to a global rotation $\phi$ around its axis. A similar conclusion holds for the third-harmonic approximation.

Finally, from Eq. (17) [or Eq. (18)] one can estimate the curvature at the loop apex

[TABLE]

and compare it with the corresponding WLC loop curvature $\kappa_{\max}=\pi c_{1}/L$ [Eq. (14)], derived in Sec. III. $\widetilde{\kappa}_{\max}$ depends on $\phi$ and, using some simple algebra, one can show that the ratio $\widetilde{\kappa}_{\max}/\kappa_{\max}$ is bounded within the interval

[TABLE]

(note that the same expression is valid both for the first- and third-harmonic approximation). Equation (21) shows that $\widetilde{\kappa}_{\max}$ does not, in general, coincide with the WLC loop apex curvature $\kappa_{\max}$. The latter has been obtained for a perfectly planar loop, however the solution (5) [and, thus, also (17) and (18)] describes an almost-planar curve with small off-planar oscillations Skoruppa et al. (2018), which are induced by the combined effect of bending anisotropy and twist-bend coupling.

IV.2 Numerical analysis

In order to test the validity of Eqs. (17) and (18), we performed Monte Carlo simulations of a “triad model”, which is derived from a direct discretization of Eq. (3). A DNA molecule with $N$ base pairs is represented by $N$ beads, each carrying a set of three orthogonal unit vectors forming the triad $\{\mathbf{\widehat{\bf e}}_{1},\mathbf{\widehat{\bf e}}_{2},\mathbf{\widehat{\bf e}}_{3}\}$, with $\mathbf{\widehat{\bf e}}_{3}$ being the tangent, and hence pointing towards the next bead. Consecutive beads are separated by a fixed distance $a=0.34$ nm, corresponding to the average base pair distance of DNA. The simulations are performed at sufficiently-low temperature, so that the system converges to its lowest-energy state (more details can be found in Ref. Caraglio et al. (2019)). Figure 2 shows the bending densities ($\Omega_{1}$, $\Omega_{2}$) and excess twist density ($\Omega_{3}$) as functions of the rescaled arc-length coordinate $s/L$ obtained from Monte Carlo simulations of a loop of $104$ bp (black circles). The lines are plots of Eqs. (17) (red) and Eqs. (18) (blue). The only adjustable parameter is the global phase $\phi$, as the stiffness constants $A_{1}$, $A_{2}$, $C$ and $G$ are input parameters (see caption of Fig. 2), while $c_{1}$ and $c_{3}$ are the universal constants given in Table 1. Figure 2 shows an excellent agreement between both harmonic approximations and the Monte Carlo data. As in the case of WLC loops shown in Fig. 1, the first-harmonic approximation overestimates the curvature at the loops ends (see inset of Fig. 2). The maximal excess twist in a loop with $104$ bp is $\max|\Omega_{3}|\approx 0.1~{}\text{nm}^{-1}$, which, compared to the average intrinsic twist $\omega_{0}=1.75~{}\text{nm}^{-1}$, corresponds to a deviation of 6% from $\omega_{0}$. Finally, note that the maximum curvature $\kappa_{\text{max}}=0.24$ nm*-1*, predicted by Eq. (14), is once more substantially lower than the total curvature in the middle of the loop of the triad model, which from the numerical data is found to be $\widetilde{\kappa}_{\max}=0.34$ nm*-1*. This is due to off-planar oscillations along the loop, as discussed above [Eq. (21)].

To further corroborate the harmonic-loop approximations, we performed low-temperature ($T=1$ K) computer simulations of 100-bp oxDNA loops (Fig. 3). OxDNA is a coarse-grained model, which describes DNA as two intertwined strands of rigid nucleotides Ouldridge et al. (2010). For the simulations we used the latest version oxDNA2 Snodin et al. (2015), which was recently found to have a substantially-nonzero twist-bend coupling constant Skoruppa et al. (2017). The electrostatic interactions are implicitly modelled through a Debye-Hückel potential. The initial configuration was a torsionally-relaxed helix, with the molecular axis having the shape of Eq. (12). To avoid undesired interactions between the two DNA ends, we constrained them at a nonzero distance $d=L/10$ by means of strong harmonic bonds. These bonds connect the center of mass of the fifth base pair at each end (yellow beads in Fig. 3b) with a fixed point in space (distance $d$ between the two centers of mass). This is to prevent end-point denaturation effects, which can affect in particular the room temperature simulations, such as those presented in Section V. The simulations were performed with the recently-developed LAMMPS Plimpton (1995) implementation of the oxDNA model Henrich et al. (2018). Figure 3a shows $\Omega_{i}$ as functions of the arc-length parameter, and once more reveals an excellent agreement with the harmonic approximations. Again, the only free parameter is the global phase $\phi$, as the stiffness constants $A_{1}$, $A_{2}$, $C$ and $G$ have already been calculated for oxDNA2 in Ref. Skoruppa et al. (2017), while $c_{1}$ and $c_{3}$ were taken from Table 1 with $d=L/10$.

V Effect of Thermal Fluctuations

Having analyzed the minimal-energy conformation of DNA loops, we now consider the effect of thermal fluctuations. For this purpose, we simulated oxDNA loops using the same setup as discussed in Section IV, with the only difference that the temperature now is $T=300$ K. Figure 4 shows kymogram traces of $\Omega_{i}$ along the DNA loop as functions of time. To better distinguish regions of predominantly-positive $\Omega_{i}$ from those of negative $\Omega_{i}$, we applied a Gaussian filter on both the spatial ($3$-bp variance) and temporal direction (4-frame variance), and used a binary color code (white for negative and blue for positive values). At the loop ends, the image is more blurred due to the thermal fluctuations dominating over the wave amplitude, as expected from the vanishing curvature. In the central region, however, clear wave patterns are visible in the bending and twist deformations, with a period following that of the DNA double helix, in line with the ground-state solutions (17) and (18). In order to better compare the relative phase among $\Omega_{i}$, the side figure shows a zoom-in of the three variables at a given small time interval (red vertical lines), further confirming the validity of Eqs. (17) and (18). We, thus, conclude that, even under the effect of thermal fluctuations, not only do the bending and twist degrees of freedom retain their modulated-wave shape, but also preserve their relative phase difference.

The global phase $\phi$ at the bottom of Fig. 4 was obtained by keeping track of the orientation of the two loop ends, which yielded the phase difference $\Delta\phi$ between successive time steps. The initial phase was determined from the Fourier analysis of the $\Omega_{2}$ wave of the first frame. The mean-squared displacement of $\phi$ grows linearly with time, as shown in Fig. 5, yielding a diffusion constant of $D=121~{}\text{rad}^{2}/\mu s$. Note that, the simulation time scale depends on the value of the Langevin damping parameter, which in this case was $6.06$ ps. The origin of the observed diffusive behavior stems from the very weak contribution of $\Delta E(\phi)$ to Eq. (19), as discussed in Section IV, indicating that the phase $\phi$ moves on an almost-flat energy landscape.

VI Conclusion

In this paper, we have investigated the shape and dynamics of short DNA loops, i.e., consisting of about 100 bp. Such loops are generated in vivo by DNA binding proteins. We have developed a two-parameter variational explicit solution for the minimal energy shape of the loop, which we refer to as harmonic loop, being described by a combination of simple trigonometric functions. This solution was found to be in excellent agreement with numerical simulations of two different coarse-grained models, and reproduces the exact WLC energy Yamakawa and Stockmayer (1972) up to four significant digits. We focused particularly on the effect of twist-bend coupling, which is an interaction arising from the DNA grooves asymmetry Marko and Siggia (1994). As recently discussed in the case of DNA minicircles and nucleosomal DNA Skoruppa et al. (2018), we found that also in DNA loops the bending deformation induces twist waves, i.e., twist oscillations following the periodicity of the helical pitch. Differently from the minicircle case, in the loops discussed here twist waves have a modulated amplitude, which is maximal at the loop apex and vanishes at the two ends. For a loop of 100 bp the maximal degree of over- and undertwisting close to the apex was estimated to be $6\%$ relative to the intrinsic double-helix twist $\omega_{0}$.

As an alternative approach, we could have used the numerically exact solution by Yamakawa and Stockmayer Yamakawa and Stockmayer (1972) to obtain the radius of curvature, $R(s)$, of the loop as a function of the curvilinear coordinate $s$. The loop shape for the model (3) is then obtained by replacing the numerical values of $R(s)$ into Eqs. (5). This approach, however, would not allow for a direct estimate of derived quantities, such as twist oscillations, curvature variation and minimal energy. The harmonic-loop approximation provides simple, yet accurate, expressions for these quantities, see e.g., Eqs. (14), (19) and (21).

Finally, we considered the loop kinetics in oxDNA simulations at room temperature. Interestingly, the bending and twist waves were not masked by the presence of thermal fluctuations, and performed a correlated motion over the whole length of the loop. This allowed us to characterize the kinetics in terms of a single parameter $\phi(t)$, describing the absolute phase of the waves. This was found to follow a simple diffusive motion, which originates from the ground-state quasi-degeneracy in $\phi$. These findings, and particularly the simplicity of our solution, may form the basis for more complex analytical calculations that involve strongly-bent DNA, such as under the action of DNA-binding proteins.

Acknowledgements.

We acknowledge financial support from the Research Funds Flanders (FWO Vlaanderen) Grant No. VITO-FWO 11.59.71.7N and FWO-SB 1SB4219N and from KU Leuven grant C12/17/006.

Appendix A Exact variational calculus

In what follows, we briefly review the derivation of the exact solution by Yamakawa and Stockmayer Yamakawa and Stockmayer (1972). Owing to the symmetry of the problem, a sufficient condition in order to ensure that the loop will close is to require that the apex of the loop has not shifted along the x-axis. Let us define $\mathbf{R}$ the vector connecting the end-point with the loop apex. We require that

[TABLE]

The determination of the minimum energy under the above constraint can be performed by introducing a Lagrange multiplier $\mu$ as follows

[TABLE]

where $\dot{\theta}\equiv\text{d}\theta/\text{d}s$. The Euler-Lagrange equation then becomes

[TABLE]

which has the following integral of motion

[TABLE]

Interpreting $\theta$ as the coordinate of a fictitious particle with mass $A$, and $s$ as the time variable, Eq. (24) can be viewed as the equation of motion for a particle in the potential $U(\theta)=\mu\cos\theta$, describing the dynamics of a pendulum under gravity, where $\mu$ plays the role of gravitational acceleration. This is the well-known Kirchhoff kinetic analogy, showing that the static conformations of elastic rods are formally equivalent to the kinetic of spinning tops. In this analogy Eq. (25) expresses the conservation of the mechanical energy. The boundary conditions (9) imply zero velocity at the begin and end point, hence $\Gamma=\mu\cos\alpha$, with $\theta(0)=\alpha$ and $\theta(L)=2\pi-\alpha$. The trajectory $\theta(s)$ is then given by Yamakawa and Stockmayer (1972)

[TABLE]

where $k^{-1}\equiv\cos(\alpha/2)$ and $F(\phi,k)$ is the incomplete Elliptic integral of the first kind

[TABLE]

The values for $\alpha$ and $\mu$ are fixed by requiring that the loop has length $L$ and that it is closed, e.g., that $\theta(s)$ satisfies (22).

Appendix B Harmonic loop ansatz

The potential energy of the pendulum has a minimum at $\theta=\pi$, corresponding to the loop apex. As a simple approximation we expand this potential around the minimum, which gives

[TABLE]

corresponding to a pendulum in the small-oscillation limit. The solution with $\theta(L/2)=\pi$ and zero velocity at the end points is then

[TABLE]

Note that the derivative of this solution is given by Eq. (10). Plugging the above equation in (22), and with a simple change of variables, one obtains

[TABLE]

where $J_{0}$ is the zeroth-order Bessel function of the first kind 111The zeroth order Besssel function of the first kind has the following integral representation

$J_{0}(x)=\frac{1}{\pi}\int_{0}^{\pi}\cos(x\cos\phi)\,d\phi,$

see Abramowitz and Stegun (1972), p. 360. This constraint, thus, fixes the constant to the smallest zero of $J_{0}$: $c_{1}=2.4048$.

Similarly, the third-order harmonic loop is obtained from the following trial function

[TABLE]

with its derivative being given by Eq. (12), and the loop-closure constraint taking the form

[TABLE]

Imposing this constraint, together with the minimization of the elastic energy (1), we obtain the values $c_{1}=2.3703$ and $c_{3}=-0.2808$ for the constants.

The above calculation can be easily generalized to open loops whose end-points are kept at some finite distance $d$. In that case, the form of the harmonic solutions (29) and (31) remains identical, but the right-hand sides of Eqs. (30) and (32) are set to the nonzero value $-d/2L$ ($x$-projection of vector pointing from the right end of the loop to the apex). As the endpoints distance increases, the harmonic ansatz becomes a more accurate approximation of the full solution (see Fig. 6). This is because, as $d$ increases, $\theta(s)$ varies within an interval getting closer to $\theta=\pi$, hence the small-angle approximation (28) becomes increasingly more accurate. This is also reflected in the improved accuracy of the data in Table 1 upon increasing $d$.

Appendix C Calculation of the energy

In this section we present some additional details over the loop-energy calculation. For simplicity, we limit the analysis to the first-harmonic approximation, as the third-harmonic case follows the same approach. The energy density, obtained from Eq. (17), becomes (to simplify the notation we drop the superscript in $\Omega_{i}^{(1)}$)

[TABLE]

To obtain the expression in the second equality we have used Eq. (17) to eliminate $\Omega_{3}$. Interestingly, this relation shows that the energy is identical to that of a pure bending deformation with a reduced stiffness $\widetilde{A}_{2}$ instead of $A_{2}$. To proceed, we use the trigonometric identities

[TABLE]

To complete the calculation of the energy, one needs to integrate in $0\leq s\leq L$. The integral of $\cos(2\pi s/L)$ vanishes; inserting Eqs. (34) and (35) in Eq. (33) and integrating, one finally gets

[TABLE]

where we have separated the dominant contribution, obtained from the integration of the constant term in the right-hand side of Eqs. (34) and (35), from the part which depends on the phase $\phi$. The former is an extensive term, giving a contribution proportional to $L$, while the integration of the $\phi$-dependent part gives a very small boundary contribution. Finally, using the definition (7) of $l_{\text{b}}$, one sees that Eq. (36) reduces to Eq. (19) of the main text.