Overtwisting induces polygonal shapes in bent DNA

TL;DR

This study combines analytical and simulation approaches to show that overtwisting in DNA minicircles causes them to adopt polygonal shapes due to twist-bend coupling, affecting their elastic properties and experimental analysis.

Contribution

It demonstrates that twist-bend coupling induces polygonal shapes in overtwisted DNA and reveals that renormalized stiffness constants govern DNA energies, impacting experimental interpretations.

Findings

Overtwisting leads to polygonal DNA shapes.

Twist-bend coupling causes periodic high and low curvature regions.

Renormalized stiffness constants, not bare ones, govern DNA energies.

Abstract

By combining analytical results and simulations of various coarse-grained models we investigate the minimal energy shape of DNA minicircles which are torsionally constrained by an imposed over or undertwist. We show that twist-bend coupling, a cross interaction term discussed in the recent DNA literature, induces minimal energy shapes with a periodic alternance of parts with high and low curvature resembling rounded polygons. We briefly discuss the possible experimental relevance of these findings. We finally show that the twist and bending energies of minicircles are governed by renormalized stiffness constants, not the bare ones. This has important consequences for the analysis of experiments involving circular DNA meant to determine DNA elastic constants.

Tables

Table 1. Table 1: Values of the stiffness coefficients (in nm) for oxDNA1 and oxDNA2 obtained in Ref. Skoruppa  et al. , 2017 . Note that while oxDNA1 has negligible twist-bend coupling G < 0.3 𝐺 0.3 G<0.3  nm, this coupling is G = 30 𝐺 30 G=30  nm in oxDNA2. The values in the table define two set of parameters of the triad model referred to as set 𝐌 1 = ( A 1 , A 2 , C , G ) = ( 84 , 29 , 118 , 0 ) superscript 𝐌 1 subscript 𝐴 1 subscript 𝐴 2 𝐶 𝐺 84 29 118 0 {\bf M}^{1}=(A_{1},A_{2},C,G)=(84,29,118,0) and set 𝐌 2 = ( A 1 , A 2 , C , G ) = ( 81 , 39 , 105 , 30 ) superscript 𝐌 2 subscript 𝐴 1 subscript 𝐴 2 𝐶 𝐺 81 39 105 30 {\bf M}^{2}=(A_{1},A_{2},C,G)=(81,39,105,30) . The last column gives A ~ 2 subscript ~ 𝐴 2 \widetilde{A}_{2} and C ~ ~ 𝐶 \widetilde{C} as obtained from Eqs. ( 4 ) and ( 5 ), respectively.

	$A_1$	$A_2$	$C$	$G$	${\tilde{A}}_2$	$\tilde{C}$
oxDNA1	84(14)	29(2)	118(1)	0.1(0.2)	29	118
oxDNA2	81(10)	39(2)	105(1)	30(1)	30	82

Equations

$$\frac{\mathrm{d}{\hat{\mathbf{e}}}_i}{\mathrm{d}s}=\left[{\Omega }_1{\hat{\mathbf{e}}}_1+{\Omega }_2{\hat{\mathbf{e}}}_2+\left({\Omega }_3+{\omega }_0\right){\hat{\mathbf{e}}}_3\right]\times {\hat{\mathbf{e}}}_i,$$

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

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

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

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

$$\beta \hat{E}$$

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

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

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

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

$$\beta E_{\text{tr}}={\int }_0^L\left[\frac{{\mu }^2{\sin }^2({\omega }_{0}s)}{2A_1}+\frac{{\mu }^2{\cos }^2({\omega }_{0}s)}{2{\tilde{A}}_2}\right]\mathrm{d}s=\frac{\pi l_b}{R},$$

$$\beta \hat{E}\approx \beta E-{\int }_0^L\left[\mu {\Omega }_1\sin (\omega s)+\mu {\Omega }_2\cos (\omega s)+\lambda {\Omega }_3\right]\mathrm{d}s,$$

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

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

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

$$\Delta \omega =\frac{\lambda A_2}{C{\tilde{A}}_2}$$

$$\beta E_{\text{tc}}$$

$$\frac{{\kappa }^2{R}^2}{4}=\frac{{\tilde{A}}_2^2}{(A_1+{\tilde{A}}_2{)}^2}+{\Gamma }^2+\frac{A_1-{\tilde{A}}_2}{A_1+{\tilde{A}}_2}{\cos }^2(\omega s)-\frac{2\Gamma A_1\cos (\omega s)}{A_1+{\tilde{A}}_2},$$

$$\Gamma \equiv \frac{G}{2A_2}R\Delta \omega \approx \frac{G}{2A_2}\Delta Lk$$

$$\frac{{\tilde{A}}_2}{A_1+{\tilde{A}}_2}{\left(1-\frac{A_1+{\tilde{A}}_2}{A_1-{\tilde{A}}_2}{\Gamma }^2\right)}^{1/2}\le \frac{\kappa R}{2}\le \frac{A_1}{A_1+{\tilde{A}}_2}+|\Gamma |$$

$$|\Gamma |\le {\Gamma }^{\ast }\equiv \frac{A_1-{\tilde{A}}_2}{A_1}$$

$$\frac{A_1}{A_1+{\tilde{A}}_2}-|\Gamma |\le \frac{\kappa R}{2}\le \frac{A_1}{A_1+{\tilde{A}}_2}+|\Gamma |$$

$${\kappa }_{//}R$$

$${\kappa }_{\perp }R$$

Full text

I Introduction

At mesoscopic length scales the elastic response of double stranded DNA to mechanical stresses is usually described by the twistable worm like chain (TWLC), which is characterized by just two elastic stiffnesses corresponding to bend and twist deformations respectively Marko (2015); Moroz and Nelson (1997). At length scales of several base pairs, relevant to protein-DNA interactions, the sequence dependent geometry of the double helix leads to a rich spectrum of elastic properties Lankaš et al. (2003); Lavery et al. (2009); Pasi et al. (2014). Even at those length scales, however, it is useful to employ simplified, analytically treatable representations that allow for the identification of generic features that are otherwise masked by sequence dependent variations. A suitable model that takes into account the general features of DNA geometry while still fulfilling the requirement of analytical treatability was put forward in the mid-90’s by Marko and Siggia Marko and Siggia (1994), who derived it from the analysis of the molecular structure of DNA. These authors showed that the asymmetry between major and minor grooves generates a coupling between bend and twist Marko and Siggia (1994), supplementing the TWLC with an additional stiffness parameter. Symmetry implies the existence of twist-bend coupling, but it does not yield any information about the magnitude of the associated coupling constant. For a long time, the effect of twist-bend coupling has been ignored, assuming that such coupling would have a minor influence on DNA elastic properties. However, a recent comparative analysis of two very similar coarse-grained DNA models one with symmetric grooves and one with asymmetric groove has shown that the twist-bend coupling constant is comparable in magnitude to the other elastic parameters describing bending and torsional stiffness Skoruppa et al. (2017). A similar conclusion was drawn from the analysis of all atom simulations of DNA Lankaš et al. (2003).

Twist-bend coupling has been shown to influence elastic properties of DNA both at long Nomidis et al. (2017, 2018) and short Skoruppa et al. (2018) molecules. For instance, in $\approx 100$ base pairs bent DNA fragment twist-bend coupling induces sinusoidal standing waves in the twist Skoruppa et al. (2018). These waves are comparable in shape and magnitude to those observed in nucleosomal DNA, which is wrapped around histone proteins Richmond and Davey (2003). In longer kilobase pair DNA subject to a stretching force and to a torque, as in single molecule experiments, twist-bend coupling has been shown to lead to a rescaling of the elastic parameters Nomidis et al. (2018). The aim of this paper is to focus on minimal energy shapes of DNA minicircles of about $10^{2}$ base pairs which are overtwisted. We show, by combining analytical and numerical results, that twist-bend coupling induces distinctive shapes of DNA in which the curvature is periodically modulated in an alternance of high and low bending regions. The periodicity is close to that of a straight double helix, but depends on the degree of overtwisting.

Thermal fluctuations strongly influence the conformation of linear DNA molecules of lengths exceeding the bending persistence length, which is approximately $50$ nm. In shorter, constrained and highly bent DNA, a situation of relevance in DNA-proteins complexes, thermal fluctuations do not influence strongly the shape of the molecule, which in first approximation assumes its minimal energy conformation. Elastic energy minimization has indeed been used to obtain the shape of DNA looping out of a Lac operon Balaeff, Mahadevan, and Schulten (1999) or for DNA wrapped around histone proteins Mohammad-Rafiee and Golestanian (2005). Although in this paper we focus on minimal energy shapes of free standing over and undertwisted minicircles, we expect that our analysis will be of particular relevance for short constrained DNA, i.e. partially bound to proteins, where thermal fluctuations have a small effect. For convenience we focus here on homogeneous minicircles, which are simpler to simulate and to describe analytically. We expect, however, that the shapes discussed here also applies to other more complex cases, such as DNA looping where translational invariance is broken.

II DNA elasticity

We briefly recall the formalism used to describe a twistable polymer such as DNA (see also e.g. Ref. Brackley, Morozov, and Marenduzzo, 2014 for more details) and review some properties of the model with twist-bend coupling, focusing on the bending and torsional persistence lengths. The conformation of an inextensible twistable elastic rod can be parametrized by three strain fields $\Omega_{1}(s)$, $\Omega_{2}(s)$ and $\Omega_{3}(s)$, where $0\leq s\leq L$ is a curvilinear coordinate and $L$ the total length of DNA. $\Omega_{1}$ and $\Omega_{2}$ are bending densities also referred to as “tilt” and “roll” deformations, while $\Omega_{3}$ is the excess twist density. DNA has an intrinsic twist density $\omega_{0}=2\pi/l_{0}\approx 1.85\leavevmode\nobreak\ \text{nm}^{-1}$, with $l_{0}=3.4\leavevmode\nobreak\ \text{nm}$ corresponding to the distance of one helical turn. Given the $\Omega_{i}(s)$ the three dimensional shape of the molecule can be obtained by solving the differential equations $i=1,2,3$:

[TABLE]

where $\{{\widehat{\bf e}}_{1}(s),{\widehat{\bf e}}_{2}(s),{\widehat{\bf e}}_{3}(s)\}$ defines an orthonormal triad (the Darboux frame) where ${\widehat{\bf e}}_{3}$ is tangent to the curve, while ${\widehat{\bf e}}_{1}$ and ${\widehat{\bf e}}_{2}$ lie on the plane of the ideal, planar Watson-Crick base pairs Marko and Siggia (1994), see Fig. 1. By convention ${\widehat{\bf e}}_{1}$ is directed along the symmetry axis of the two grooves and points in the direction of the major groove. Finally ${\widehat{\bf e}}_{2}={\widehat{\bf e}}_{3}\times{\widehat{\bf e}}_{1}$, which yields to a vector connecting the two DNA backbones Olson et al. (2001). The conformation with all $\Omega_{i}=0$ describes a straight twisted rod with an intrinsic twist density $\omega_{0}$.

From the analysis of the molecular symmetry of DNA, Marko and Siggia Marko and Siggia (1994) derived the following free energy functional to lowest order in $\Omega_{i}$:

[TABLE]

where $\beta=1/k_{B}T$ is the inverse temperature, and $A_{1}$, $A_{2}$, $C$ and $G$ are the stiffness parameters. The last term in (2) couples twist ($\Omega_{3}$) with the bending towards the DNA grooves ($\Omega_{2}$). In the case of vanishing twist-bend coupling constant $G=0$ one recovers the twistable wormlike chain (TWLC) which is the commonly employed model to describe DNA elasticity. In this paper we focus in particular on effects of twist-bend coupling on the minimal energy shapes of bent DNA. Before entering into the details we recall in this Section some properties of the model (2), which differ from those of the standard TWLC.

The bending persistence length $l_{b}$ of the model (2) has been computed from the decay of the tangent-tangent correlation and it is given by the following relation Nomidis et al. (2017); Skoruppa et al. (2017)

[TABLE]

where

[TABLE]

can be viewed as a rescaled bending stiffness along the "easy" bending axis (recall that $A_{2}<A_{1}$). In the isotropic TWLC limit $A_{1}=A_{2}=A$ and $G=0$ Eq. (3) reduces to the well-know results $l_{b}=A$, i.e. $A$ is both the bending stiffness and the bending persistence length. In the anisotropic case $A_{1}\neq A_{2}$ and $G=0$, Eq. (3) shows that $l_{b}$ is the harmonic mean of $A_{1}$ and $A_{2}$, a known result, see Ref.s Lankaš et al., 2000; Eslami-Mossallam and Ejtehadi, 2008. In the general case $G\neq 0$ the result can be cast in the form of a harmonic mean of $A_{1}$ and $\widetilde{A}_{2}$. From the analysis of the twist correlations one can calculate the twist correlation length Nomidis et al. (2017) and finds $l_{t}=2\widetilde{C}$ with

[TABLE]

which reduces to the well-known result $l_{t}=2C$ in the TWLC limit $G=0$. Equations (3) and (5) are exact and have been obtained from the analysis of thermal fluctuations of the model (2). These quantities are also relevant in minimal energy DNA conformations that are permanently bent and twisted, as we will show in the next Sections.

III DNA minicircles

A considerable amount of studies have been devoted to the analysis of equilibrium and kinetic properties of DNA minicircles Mohammad-Rafiee and Golestanian (2003); Olson, Swigon, and Coleman (2004); Du et al. (2005); Mohammad-Rafiee and Golestanian (2005); Lankaš, Lavery, and Maddocks (2006); Norouzi, Mohammad-Rafiee, and Golestanian (2008); Demurtas et al. (2009); Lionberger et al. (2011); Wang and Pettitt (2016); Schurr (2017); Pasi et al. (2017) or to short DNA loops obtained by bending two ends of DNA Sankararaman and Marko (2005); Mulligan et al. (2015). Our aim is to describe the effect that twist-bend coupling has on the minimal energy shape of overtwisted minicircles. We employ a homogeneous model (2), neglecting sequence-dependent variations of the elasticity, which we expect to be valid when averaging over different sequences. In a minicircle, the two ends of the double helix are covalently sealed and in order for the loop to close one usually needs to over or underwind the molecule so that the end point meet in phase. An excess linking number can also be induced by using appropriate enzymes that overtwist the molecule. In what follows we discuss separately the torsionally relaxed and the over and under-twisted cases.

III.1 Torsionally relaxed minicircles

As this case has been discussed recently in Ref. Skoruppa et al., 2018, we just briefly review it here. An analytical shape for anisotropic bending stiffness and $G\neq 0$ cannot be found easily. The constraint to be imposed to form circular DNA requires that the end points of the strands meet each other smoothly and that the tangent vector is continuous. This constraint is typically expressed using lab-frame quantities, as Euler angles describing the configuration of the DNA with respect to some fixed axes. The minimal energy shapes could then be calculated numerically. Here, however, we employ an approximation that allows us to perform a local energy minimization, expressing the constraint in terms of $\Omega$’s. We assume that the exact solution is a small perturbation of a perfect circle and that the solution is periodic over helical repeats.

The vector $\vec{\Omega}_{b}\equiv\Omega_{1}{\widehat{\bf e}}_{1}+\Omega_{2}{\widehat{\bf e}}_{2}$ is the bending density. Given a fixed unit vector ${\widehat{\bf x}}$ (see Fig. 1), Ref. Skoruppa et al., 2018 introduced the following local constraint $-\mu\leavevmode\nobreak\ \vec{\Omega}_{b}\cdot{\widehat{\bf x}}$, with $\mu$ a Lagrange multiplier. This term favors the alignment of $\vec{\Omega}_{b}$ along ${\widehat{\bf x}}$, i.e. it forces the molecule to be bent and to remain close to a plane orthogonal to ${\widehat{\bf x}}$ (see Fig. 1). For a straight double helix in the plane orthogonal to ${\widehat{\bf x}}$ and with constant twist density equal to $\omega_{0}$, one can choose the curvilinear coordinate such that ${\widehat{\bf x}}=\sin(\omega_{0}s){\widehat{\bf e}}_{1}+\cos(\omega_{0}s){\widehat{\bf e}}_{2}$. This relation remains approximately valid if within one helical turn bending is weak and local twist variations are small, conditions which can be respectively expressed as $|\vec{\Omega}_{b}|\ll\omega_{0}$ and $|\Omega_{3}|\ll\omega_{0}$. Within this approximation the constraint takes the following form:

[TABLE]

with $E$ given by (2). Minimization with respect to $\Omega_{i}$ gives the following solution for a torsionally relaxed minicircle Skoruppa et al. (2018):

[TABLE]

In the previous equations ${\mu}\equiv l_{b}/R$, where $R$ is the average radius of the circle and where $l_{b}$, given by (3), is the bending persistence length of the model (see Ref. Skoruppa et al., 2018 for more details).

Equations (7) yield a perfect circle only in the case of $A_{1}=\widetilde{A}_{2}$, corresponding to a constant curvature $\kappa=\sqrt{\Omega_{1}^{2}+\Omega_{2}^{2}}$. In all other cases the above equations describe a quasi-circular shape that exhibits small off-planar and in-planar oscillations, with a modulated total curvature $\kappa$. The analysis of Eq. (7) shows that the dimensionless quantity $\kappa R$ ($R$ is the average radius of curvature of the minicircle) is bounded within the interval

[TABLE]

which reduces to a perfect circle $\kappa R=1$ in the isotropic TWLC limit $A_{1}=A_{2}$, $G=0$ (so that $\widetilde{A}_{2}=A_{2}$). Note that $A_{2}$ is the bending stiffness along the easy axis hence $A_{2}<A_{1}$ and twist-bend coupling enhances the curvature anisotropy since ${\widetilde{A}}_{2}<A_{2}$. Twist-bend coupling induces oscillations in the twist, referred to as "twist waves" in Ref. Skoruppa et al., 2018, which are in antiphase with the "roll" ($\Omega_{2}^{(tr)}$) wave. Twist oscillations are experimentally observed in crystal structures of DNA wrapped around histone proteins Richmond and Davey (2003), however their origin has been so far attributed to an interaction with the underlying histone core proteins, while Eq. (7) shows that these oscillations are a natural feature of bent DNA, directly deriving from the effect of twist-bend coupling Skoruppa et al. (2018).

The elastic energy associated with this configuration is obtained by inserting equations (7) into (2). Using $\mu=l_{b}/R$ and $L=2\pi R$, for a torsionally relaxed minicircle one finds

[TABLE]

which is formally identical to the energy of a TWLC minicircle. The difference being that in the TWLC the persistence length $l_{b}$ is the harmonic mean of the bending stiffnesses, whereas in model (2) the same quantity is a function of all elastic parameters of the model, see (3) and (4).

III.2 Over and undertwisted minicircles

In the case of torsionally constrained minicircles, one should impose a constraint of a fixed linking number $Lk$. The White theorem states that the linking number is the sum of twist and writhe $Lk=Tw+Wr$. As we are interested in quasi-planar conformations the writhe is small and $Lk\approx Tw$. Therefore, we impose a constraint on a twist instead by introducing an additional Lagrange multiplier as follows

[TABLE]

where we have used the same approximation as in (6) but with ${\widehat{\bf x}}=\sin(\omega s){\widehat{\bf e}}_{1}+\cos(\omega s){\widehat{\bf e}}_{2}$, in order to take into account that the introduction of a constraint in $\Omega_{3}$ induces a shift in the intrinsic twist to $\omega\equiv\omega_{0}+\Delta\omega$. Local energy minimization yields the modified equations

[TABLE]

As expected, a non-zero $\lambda$ introduces an offset in $\Omega_{3}$, i.e. an average excess twist density given by

[TABLE]

where again we assumed that the writhe has negligible contribution. Because of twist-bend coupling there is a corresponding offset in $\Omega_{2}$ as well, which enhances the inhomogeneity in the curvature. The shape described by Eqs. (11) alternates between high and low curvature regions, but the excess twist induces oscillations in $\kappa R$ in a wider range when compared to the torsionally relaxed case (8) (more details are shown in Appendix A). The $\Omega$’s from Eqs. (11) are in very good agreement with numerical results on coarse-grained DNA models discussed in the next section.

Similarly to what was done for the untwisted case Skoruppa et al. (2018), one can show that the Lagrange multiplier $\mu$ fixes the average radius of curvature, hence $\mu=l_{b}/R$. Plugging Eqs. (11) into (2) we obtain the following energy for a torsionally contrained minicircle

[TABLE]

where $\widetilde{C}$ has been defined by Eq. (5) and the term proportional to $\cos(\omega s)$ averages out in the integration. As in the untwisted case the energy is formally identical to that of a TWLC. Again in model (2), the difference with TWLC is that the bending stiffness $l_{b}$ and torsional stiffness $\widetilde{C}$ are function of all elastic constant. Naturally, for large excess twist the approximation of negligible writhe will break down and the minicircle will start supercoiling. This instability is not captured by the present model (10), as it was derived by a local constraint, explicitly neglecting this effect.

IV Numerical Results

In order to check the theoretical predictions, we performed some numerical calculations using two different models. The first one, which we refer to as the triad model (see also Refs. Skoruppa et al., 2018; Nomidis et al., 2018) is obtained by a direct discretization of the continuum model (2). The second model is oxDNA Ouldridge, Louis, and Doye (2010), a coarse-grained model describing DNA as two intertwined strings of rigid nucleotides. In the triad model the stiffness constants $A_{1}$, $A_{2}$, $C$ and $G$ are input parameters and can be freely chosen, while for oxDNA their values are fixed by the force field parametrization, which was tuned to reproduce known structural, thermodynamical and mechanical properties of DNA Ouldridge, Louis, and Doye (2010). In addition, oxDNA comes in two versions: an implementation with symmetric grooves (oxDNA1) and an improved version that explicitly introduces asymmetric grooves (oxDNA2). This feature makes the model suitable to test the implications of twist-bend coupling, which derives precisely from the groove asymmetry Marko and Siggia (1994). The elastic parameters for oxDNA1 and oxDNA2 were computed in Ref.Skoruppa et al., 2017 from the analysis of equilibrium fluctuations of a linear molecule and the results are shown in Table 1, yielding in particular a value of $G$ which is comparable to that of the other elastic constants.

IV.1 MC simulations with the Triad Model

In the triad model, a double stranded DNA of $N$ base pairs is represented by $N$ beads, each carrying a frame of three orthogonal unit vectors $\{\widehat{\bf e}_{1},\widehat{\bf e}_{2},\widehat{\bf e}_{3}\}$. The distance between consecutive beads is fixed and equal to $a=0.33$ nm and the vector $\widehat{\bf e}_{3}$ always points towards the sequentially adjacent bead. Given two consecutive triads, the deformation parameters are defined by a definition analogous to Eq. (1) valid for finite rigid body rotations (for details see Supplemental Material of Ref. Skoruppa et al., 2017), while the energy of a conformation is obtained by discretizing the continuum energy model, Eq. (2). Low temperature Monte Carlo (MC) simulations have been carried out using the two sets of parameters ${\bf M}^{1}$ and ${\bf M}^{2}$ matching the oxDNA1 and oxDNA2 values (see caption of Table 1).

Figure 2 shows a comparison of the $\Omega_{i}$ obtained from the triad model simulations (colored lines) with the analytical predictions (black lines). In the simulations we fix the linking number $Lk$, which is a topological invariant measuring the number of times the two strands are wound around each other. The White theorem states that the linking number is the sum of twist and writhe $Lk=Tw+Wr$. As we are interested in quasi-planar conformations the writhe is small and $Lk\approx Tw$. Therefore the constraint used in (10) to fix the excess twist density $\Omega_{3}$ is expected to adequately describe closed circular DNA. The left panel of Fig. 2 shows a torsionally relaxed minicircle with $104$ base pairs and linking number $Lk=10$, while in the right panel the circle is overtwisted and has $Lk=11$. All analyzed cases display excellent agreement between analytical models and MC data. The horizontal axis shows a single helical turn, corresponding to a phase $0\leq\phi\leq 2\pi$ obtained from the analysis of the Fourier spectrum of $\Omega_{i}$, as explained in the caption. In a torsionally relaxed DNA one helical turn corresponds to $10.4$ base pairs as shown in the horizontal top scale of the left panel. In the overtwisted case ($Lk=11$) one helical turn corresponds to $9.5$ base pairs.

As shown in Fig. 2, in the torsionally relaxed case $Lk=10$ and for the set ${\bf M}^{2}$ ($G\neq 0$) the twist oscillates, while these oscillations are absent for the set ${\bf M}^{1}$ ($G=0$). In the torsionally constrained case $Lk=11$ the twist of both sets is shifted, see bottom right graph of Fig. 2. Overtwisting affects the $\Omega_{2}$ of set ${\bf M}^{2}$, but not that of set ${\bf M}^{1}$, as predicted by the analytical model. Is is worth emphasizing that the theoretical predictions of both panels do not contain adjustable parameters. In fact, the Lagrange multipliers $\mu$ and $\lambda$ are fixed as follows: $\mu=l_{b}/R$, where $R$ is the radius of a perfect circular chain of $104$ beads with a distance between consecutive beads of $0.33$ nm ($R\simeq 5.46$ nm), $\lambda=C\tilde{A}_{2}\Delta\omega/A_{2}$ with $\Delta\omega=0.18$ nm*-1*, and the stiffness parameters are given. In the undertwisted case (not shown) very similar profiles for the $\Omega_{i}$ are observed, but the shifts in $\Omega_{2}$ and $\Omega_{3}$ carry the opposite sign.

Figure 3 shows the typical shapes of relaxed and overtwisted minicircles. In both sets ${\bf M}^{1}$ and ${\bf M}^{2}$, the relaxed configuration ($Lk=10$) is almost completely planar and closely resembles a perfect circle. However, the introduction of additional twist leads to remarkable differences in the behavior of minicircles parameterized by the sets ${\bf M}^{1}$ and ${\bf M}^{2}$. The latter starts to exhibit the shape of a rounded polygon or, more precisely, a rounded hendecagon, where the amount of vertices is induced by the imposed excess linking number. Furthermore, the structure is moderately off-planar, as illustrated by the plot of $\tilde{d}$, the signed distance of each base pair from the best fitting plane vs. base pair position not (see left lower panel of Fig. 3). Consistent with the oscillations of the strain fields $\Omega_{i}$, $\tilde{d}$ fluctuates with a wavelength of $2\pi/\omega$. On the other hand, when one considers the parameter set without twist-bend coupling (${\bf M}^{1}$), depending on the setup of the MC simulation, two typical configurations are found shown in blue and green in Fig. 3, right. If one starts from a perfectly circular overtwisted shape and performs MC updates at low temperature, the simulation relaxes to the shape drawn in blue (Fig. 3, right). From this shape the $\Omega_{i}$ shown in Fig. 2 were calculated. If, on the other side, one starts from a high temperature simulation and gradually lowers the temperature to reach the ground state a strongly off planar conformation, as that shown in green in Fig. 2, is obtained. The latter shows an onset of supercoiling, which is not found in the simulations with set ${\bf M}^{2}$ for the same value of $Lk$. In that case the polygon shape is always recovered at low temperatures, irrespectively from the simulation path followed. This shows that the model ${\bf M}^{1}$ is more prone to supercoiling compared to model ${\bf M}^{2}$, for which the torsional stress is released in bending and off-planar fluctuations. This is also reflected in Eq. (LABEL:eq:en_tc), which shows that the torsional energy is controlled by the parameter $\widetilde{C}$ (5), rather than the intrinsic twist stiffness $C$. The sets ${\bf M}^{1}$ and ${\bf M}^{2}$ have comparable twist stiffness $C$, but $\widetilde{C}=C=118$ nm for set ${\bf M}^{1}$ as $G=0$, while $\widetilde{C}=82$ nm and $C=105$ nm for set ${\bf M}^{2}$ which implies a considerable lower twisting energy for the same amount of overtwisting. This explains why model ${\bf M}^{1}$ supercoils more easily when compared to ${\bf M}^{2}$. Also the model ${\bf M}^{2}$ eventually supercoils at larger $|\Delta Lk|$ (not shown). We will not discuss the properties of the supercoiling transition for the two models here, which would require sampling both systems at experimentally relevant temperatures. Appendix B provides some details of MC simulations with triad model, confirming that model ${\bf M}^{2}$ has a lower propensity to supercoiling compared to ${\bf M}^{1}$ once the same external parameters as linking number and temperatures are chosen.

IV.2 MD simulations with oxDNA

Double helical coarse-grained models have become very popular in the recent few years to study a large number of equilibrium and dynamical properties of DNA Ouldridge, Louis, and Doye (2010); Šulc et al. (2012); Fosado et al. (2016); Frederickx, In’t Veld, and Carlon (2014); Lequieu, Schwartz, and de Pablo (2017); Li and Kabakçıoğlu (2018); Chakraborty, Hori, and Thirumalai (2018); Coronel, Suma, and Micheletti (2018). oxDNA Ouldridge, Louis, and Doye (2010) provides an effective mesoscopic description of DNA in which each nucleotide is represented by a rigid object equipped with three interactions sites for base pairing, coaxial stacking, electrostatic and steric interactions, that was tuned to reproduce the properties of dsDNA. Langevin dynamics of the system was integrated at low temperature ($15$ K) with the LAMMPS package Plimpton (1995) using the implementation of Henrich et al. Henrich et al. (2018) and default values for the interaction parameters.

Again, the $\Omega_{i}$ calculated from oxDNA MD simulations are consistent with those predicted by the analytical model (see Fig. 4). As already pointed out in Ref. Skoruppa et al., 2018, small deviations can be noticed in the amplitude of $\Omega_{1}$ oscillations. Interesting enough, also the value of $\Omega_{3}$ for oxDNA2 slightly deviates from the analytical curve but only in the case in which the minicircle is overtwisted. As for $\Omega_{1}$, these small deviations probably arise from some additional interactions (e.g. higher-order terms) present in oxDNA, but not considered in the energy functional. Notice that we considered an overtwisted minicircle with small additional twist $\Delta\omega/\omega_{0}\simeq 5$%. Such choice is due to the fact that for higher additional twist, e.g. of the order of $10$% as imposed in the case of Fig. 2, both oxDNA1 and oxDNA2 show a behavior of strong deviations from planarity and a propensity to form supercoils.

Next, we considered strongly overtwisted minicircles in confined geometries so that the DNA cannot easily supercoil. Two cases were analyzed: a DNA minicircle confined between two flat surfaces and a minicircle wrapped around a cylinder, shown in Figure 5. These two situations can be relevant, respectively, in the case of AFM experiments where DNA is confined in 2D by absorbption on a mica surface Rivetti, Guthold, and Bustamante (1996); Valle et al. (2005); Wiggins et al. (2006); Vanderlinden and De Feyter (2013) or in the case of nucleosomes where DNA is wrapped around the cylindrically shaped octamer of histone proteins Eslami-Mossallam, Schiessel, and van Noort (2016); Richmond and Davey (2003). Figure 5 shows the typical shapes of confined oxDNA1 and oxDNA2 minicircles when an excess linking number of $20$% is introduced. Apart obviously from the twist-bend coupling $G$, oxDNA1 and oxDNA2 have similar elastic parameters (see Table 1) but their response to overtwisting and confinement is very different. This is in agreement to what was found in the triad model: as for the set ${\bf M}_{1}$, also oxDNA1 has a stronger tendency to develop off-planar conformation, indicating that it is easier to supercoil. Furthermore, in the case of oxDNA2, a closer inspection of the base pairs’ center of mass, as shown in the inset of Fig. 5 (left), indicates that the minicircle tends to have a rounded polygonal shape, with the periodicity of the double helix, in agreement with the theory discussed in this paper. On the other hand, oxDNA1 has a smoother curvature even if, under planar confinement, the shape of oxDNA1, once projected in 2D, resembles a rounded square. The origin of these seemingly strong curved regions is the onset of a buckling transition, i.e. the apparent strong curvature in the four corners of the 2D projection is due to writhe rather than curvature and does not exhibit the periodicity of the inert double helical repeat length.

V Conclusion

In this paper we have studied minimal energy shapes of torsionally constrained circular DNA molecules. As shown earlier Skoruppa et al. (2018), the effect of twist-bend coupling is to produce shapes characterized by coupled oscillations in twist ($\Omega_{3}$) and in the bend ($\Omega_{2}$) densities. We have extended here the investigation to the effect of a torsional strain which induces a net shift in the twist density $\Omega_{3}$, forcing it to oscillate around a non-zero average value. As a consequence of twist-bend coupling this effect is transmitted to the groove-bending strain $\Omega_{2}$. The breaking of the symmetry of $\Omega_{2}$ oscillations results in a shape resembling that of a rounded regular polygon with a periodic alternation of high and low curvature regions. We have shown (extending the theory of Ref. Skoruppa et al., 2018) that a simple analytical model reproduces very well the shapes obtained from simulations. The comparison between theory and simulations is remarkable as there are no adjustable parameters.

The analytical model provides a simple way to estimate the energy of torsionally relaxed and torsionally constrained minicircles. It turns out that, due to the peculiar shapes of the circles induced by twist-bend coupling, the elastic energy due to bending and twist is not governed by the intrinsic stiffnesses, but by rescaled parameters $l_{b}$ and $\widetilde{C}$, given by Eq. (3) and (5) respectively. As a consequence, although the two sets of parameters used in this work have comparable torsional stiffness $C$ (see Table. 1), their torsional response is very different. In the set ${\bf M}^{1}$, with vanishing twist bend coupling, the twist energy is $\frac{1}{2}CL\Delta\omega^{2}$ while in ${\bf M}^{2}$ this is reduced to $\frac{1}{2}\widetilde{C}L\Delta\omega^{2}$ (recall that $\widetilde{C}<C$, see (5)). In the set ${\bf M}^{2}$ the minimal energy shape exploits the presence of a cross term $G\Omega_{2}\Omega_{3}$ which can become negative, hence lowering the energy, if $\Omega_{2}$ and $\Omega_{3}$ have opposite signs. The value of torsional elastic constant for DNA has been discussed at length in the literature with different techniques Horowitz and Wang (1984); Fujimoto and Schurr (1990); Vologodskii and Marko (1997); Moroz and Nelson (1998); Bouchiat and Mézard (1998); Kauert et al. (2011); Bryant, Oberstrass, and Basu (2012) providing values ranging typically from $C=75$ nm to $C=110$ nm, although occasionally lower or higher values have been reported (a table collecting the elastic constant measurements reported in the literature from various experimental techniques can be found in the supplemental of Ref. Nomidis et al., 2017). One way of extracting $C$ is from the analysis of dynamical or equilibrium properties of DNA minicircles, see eg. Refs. Shore and Baldwin, 1983; Horowitz and Wang, 1984. In this analysis it is assumed that DNA is described by an elastic rod model with independent twist and bending deformations (TWLC) and that if the circles are sufficiently small ($\lesssim 200$ bp) the contribution from writhe fluctuations can be neglected Horowitz and Wang (1984). Measurements of topoisomers, i.e. sequences of equal length that differ only by their linking number, have been used to estimate the bare torsional stiffness $C$. We have shown here that, in presence of twist-bend coupling, the energetic behavior of minicircles is governed by a renormalized stiffness $\widetilde{C}$, which is smaller than $C$ and contains the parameters $G$ and $A_{2}$, see (5).

In a recent paper Skoruppa et al. (2018) some of us showed that minimal energy shapes of minicircles with twist-bend coupling fit well structural nucleosomal DNA data, as obtained from X-ray crystallography. Nucleosomal DNA is wrapped around the nucleosome, a stable complex formed by histone proteins tightly bound to each other. The nucleosome is known to slide along the DNA and one of the most discussed mechanisms is that of the diffusion of twist defects, see Ref. Eslami-Mossallam, Schiessel, and van Noort, 2016 for a recent review. We have shown here that local under or overtwisting of bent DNA is accompanied by a change in shape with a modulation of the local curvature, which may influence the way the twist defects propagate along the nucleosomal DNA sequence. This is an interesting issue to be considered in future work.

Acknowledgements.

Discussions with M. Laleman, J. Marko, S. Nomidis and J.M. Schurr are gratefully acknowledged. MC and ES aknowledge financial support from KU Leuven grant C12/17/006.

Appendix A Curvature

From Eqs. (11) and using $\mu=l_{b}/R$ one finds

[TABLE]

where we defined

[TABLE]

($\Gamma$ is a rescaled dimensionless excess twist density, or equivalently the excess linking number if one neglects the contribution of the writhe $2\pi R\Delta\omega\approx\Delta Lk$). The torsionally relaxed case corresponds to $\Gamma=\Delta\omega=0$, $\omega=\omega_{0}$. Equation (14) gives in this case a maximal curvature when $\cos(\omega_{0}s)=\pm 1$ and a minimal curvature when $\cos(\omega_{0}s)=0$. These are the bounds given in Eq. (8).

In the torsionally constrained case and for non vanishing twist-bend coupling one has $\Gamma\neq 0$. The analysis of (14) yields the following bounds

[TABLE]

valid for

[TABLE]

and

[TABLE]

for $|\Gamma|\geq\Gamma^{*}$. In the limit $\Gamma=0$, Eq. (16) reduces to (8). Comparing (16) and (18) with (8) one sees that introducing an excess twist indeed increases the range of values through which the curvature $\kappa$ oscillates.

To gain some more insigths on the different contributions to $\kappa$ one can decompose it into an in-planar component $\kappa_{\mathbin{\!/\mkern-4.0mu/\!}}=\vec{\Omega}_{b}\cdot{\widehat{\bf x}}$ and an off-planar component $\kappa_{\perp}=\vec{\Omega}_{b}\cdot[{\widehat{\bf e}}_{3}\times{\widehat{\bf x}}]$, where $\widehat{\bf x}$ identifies the plane where the minicircle lies (see Fig. 1). Obviously $\kappa^{2}=\kappa_{\mathbin{\!/\mkern-4.0mu/\!}}^{2}+\kappa_{\perp}^{2}$. We recall the definition of bending vector $\vec{\Omega}_{b}=\Omega_{1}\widehat{\bf e}_{1}+\Omega_{2}\widehat{\bf e}_{2}$. Using the approximation of the main text $\widehat{\bf x}=\sin(\omega s)\,{\widehat{\bf e}}_{1}+\cos(\omega s)\,{\widehat{\bf e}}_{2}$ one finds:

[TABLE]

The top panels of Fig. 6 show the total curvature $\kappa$ and in-plane curvature $\kappa_{\mathbin{\!/\mkern-4.0mu/\!}}$ as given by (14) and (19), respectively. There is a small difference between the two, showing that the contribution of off-planar bending to the total curvature, given by $\kappa_{\perp}$, is small. For an isotropic model ($A_{1}=\widetilde{A}_{2}$) with $\Gamma=0$ (corresponding to either $\Delta Lk=0$ or $G=0$) one recovers from (14), (19) and (20) a perfect homogeneous and planar circle of radius $R$: $\kappa=\kappa_{\mathbin{\!/\mkern-4.0mu/\!}}=1/R$ and $\kappa_{\perp}=0$. The introduction of anisotropic bending leads to curvature oscillations, due to the terms proportional to $\cos(2\omega s)$ and $\sin(2\omega s)$ in (19) and (20). The period of these oscillations is half a helical repeat. Overtwisting or undertwisting minicircles with non-zero twist-bend coupling ($\Gamma\neq 0$) breaks this symmetry by introducing terms that oscillate with a period of a full helical repeat length. The difference between undertwisting ($\Gamma<0$) and overtwisting ($\Gamma>0$) is that in the former case the region of maximal curvatures $\kappa$ or $\kappa_{\mathbin{\!/\mkern-4.0mu/\!}}$ corresponds to a global maximum of $\Omega_{2}$ ($>0$) while in the latter to a global minimum, see lower panels of Fig. 6. The different signs of $\Omega_{2}$ correspond to different mode of groove bendings, with a positive value of $\Omega_{2}$ corresponding to a bending towards the major groove.

Appendix B Writhe behavior in MC simulations with triad model

Figure 7 shows plots of the writhe, $Wr$, as a function of the MC time steps for the triad model. The data are for overtwisted minicircles with $N=104$ bp and $Lk=11$, obtained using the same parametrizations considered in Fig. 3 and complement the results presented that figure. Two different temperature runs are shown $T=10^{-2}$ (Fig. 7(a)) and $T=1$ (Fig. 7(b)). In this scale $T=1$ corresponds to room temperature. The simulations are performed for both sets ${\bf M}^{1}$ and ${\bf M}^{2}$ with two different initial corresponding to a perfect circular shape with $Wr=0$ and a supercoiled conformation with $Wr\simeq 1$. At $T=10^{-2}$ and in the set ${\bf M}^{1}$ the dynamics does not change significantly the writhe, which slightly increase for the circular initial condition to $Wr\simeq 0.02$. This state corresponds to the circular shape shown in Fig. 3, right (blue circle $Lk=11$, set ${\bf M}^{1}$). The writhe is roughly constant also when starting from the supercoiled initial condition corresponding to the strongly off-planar shape of Fig. 3, right (green $Lk=11$, set ${\bf M}^{1}$). The low $T$ simulations for ${\bf M}^{1}$ hence show that there are two local free energy minima corresponding to the planar and supercoiled state. For the set ${\bf M}^{2}$ and $T=10^{-2}$ both initial conditions converge to $Wr\simeq 0.05$, corresponding to the polygonal shape discussed earlier Fig. 3, left. Note that the simulation remains for long time in the supercoiled state before relaxing to the polygonal shape, indicating that a supercoiled state is metastable. The writhe is higher than that of the circle of ${\bf M}^{1}$, in accordance with the out of plane oscillations presented in Fig. 3. The analysis was extended to runs at room temperature $T=1$ (Fig. 7(b)). In this case fluctuations in the shapes are higher. For the set ${\bf M}^{1}$ both initial conditions converge to the supercoiled state, whereas for ${\bf M}^{2}$ a low writhe conformation is reached. Overall, the results confirm the low propensity towards supercoiling in the case with a non-vanishing twist-bend coupling.