Strongly bent double-stranded DNA: reconciling theory and experiment

TL;DR

This paper introduces a unified framework for understanding polymer bending, especially in DNA, that accounts for both weak and strong bending regimes and aligns well with experimental data.

Contribution

It presents a new quantitative model for polymer bending that captures the transition from harmonic to linear energy behavior, explaining DNA cyclization experiments.

Findings

The effective bending energy of DNA becomes non-convex beyond 2 degrees per base-pair.

The model predicts increased cyclization probability for very short DNA loops.

The non-convexity is explained by short-range Lennard-Jones interactions in DNA.

Abstract

The strong bending of polymers is poorly understood. We propose a general quantitative framework of polymer bending that includes both the weak and strong bending regimes on the same footing, based on a single general physical principle. As the bending deformation increases beyond a certain (polymer-specific) point, the change in the convexity properties of the effective bending energy of the polymer makes the harmonic deformation energetically unfavorable: in this strong bending regime the energy of the polymer varies linearly with the average bending angle as the system follows the convex hull of the deformation energy function. For double-stranded DNA, the effective bending deformation energy becomes non-convex for bends greater than about 2 degrees per base-pair, equivalent to the curvature of a closed circular loop of about 160 base pairs. A simple equation is derived for the polymer loop energy that covers both the weak and strong bending regimes. The theory shows quantitative agreement with recent DNA cyclization experiments on short DNA fragments, while maintaining the expected agreement with experiment in the weak bending regime. Counter-intuitively, cyclization probability (j-factor) of very short DNA loops is predicted to increase with decreasing loop length; the j-factor reaches its minimum for loops of 45 base pairs. Atomistic simulations reveal that the attractive component of the short-range Lennard-Jones interaction between the backbone atoms can explain the underlying non-convexity of the DNA effective bending energy, leading to the linear bending regime. Applicability of the theory to protein-DNA complexes, including the nucleosome, is discussed.

Tables

Table 1. Table 1 : j-factor ratios, J ​ ( L 1 ) / J ​ ( L 2 ) 𝐽 subscript 𝐿 1 𝐽 subscript 𝐿 2 J(L_{1})/J(L_{2}) , predicted using ECH and WLC models compared with experiment Vafabakhsh and Ha ( 2012 ) .

$L_1$(bp)	$L_2$(bp)	Experiment	ECH	WLC
40	50	1.50	0.993	1.12$\times {10}^{-6}$
71	101	1.51$\times {10}^{-1}$	2.01$\times {10}^{-1}$	3.08$\times {10}^{-6}$
80	101	2.17$\times {10}^{-1}$	3.28$\times {10}^{-1}$	6.10$\times {10}^{-4}$
90	101	3.56$\times {10}^{-1}$	5.59$\times {10}^{-1}$	3.73$\times {10}^{-2}$

Equations

$$E_{chain}=\sum _i^{N-1}\frac{1}{2}{k}_{B}T\frac{L_p}{l}{\theta }_i^2$$

$$E_{chain}=N\text{mathbfcal}E(\bar{\theta })=\underset{{\sum }_{i=1}^N{\theta }_i=N\bar{\theta }=\alpha }{\min }\{\sum _i^{N}E({\theta }_i)\}$$

$$\text{mathbfcal}E(\bar{\theta })=\frac{\bar{\theta }-{\theta }_a}{{\theta }_b-{\theta }_a}(E({\theta }_b)-E({\theta }_a))+E({\theta }_a)$$

$$\text{mathbfcal}E(\bar{\theta })=\{\begin{array}{ll}\frac{1}{2}{k}_{B}T{L}_p{\bar{\theta }}^2& \text{if }\bar{\theta }\le {\theta }_a\\ & \\ k_{B}T{L}_p{\theta }_a(\bar{\theta }-\frac{1}{2}{\theta }_a)& \text{if }{\theta }_a<\bar{\theta }<{\theta }_b\end{array}$$

$$\mathbfcal E_{loop}(L)=\begin{cases}2\pi^{2}k_{B}T\frac{L_{p}}{L}&\text{if }L>\frac{2\pi}{\theta_{a}}\\ \\ k_{B}TL_{p}\theta_{a}\Big{(}2\pi-\frac{1}{2}L\theta_{a}\Big{)}&\text{if }\frac{2\pi}{\theta_{b}}<L<\frac{2\pi}{\theta_{a}}\end{cases}$$

$$j(L)\simeq \frac{k}{L_p^3}{\left(\frac{L_p}{L}\right)}^5\exp \left(-\frac{\text{mathbfcal}{E}_{loop}}{k_{B}T}+\frac{L}{4L_p}\right)$$

$$j(L)\simeq \frac{k}{L_p^3}{\left(\frac{L_p}{L}\right)}^5\exp \left(-\frac{\text{mathbfcal}{E}_{loop}(L)}{k_{B}T}+\frac{L}{4L_p}\right)$$

$$J(L)$$

$$\displaystyle+J_{bot}\Big{(}1-\cos\Big{(}\frac{2\pi L}{10}\Big{)}\Big{)}\Big{]}$$

$$J_{top/bot}(L)=\frac{k_{top/bot}}{L_p^3}{\left(\frac{L_p}{L}\right)}^5{e}^{\left(L_p{\theta }_a\left(2\pi -\frac{1}{2}L{\theta }_a\right)+\frac{L}{4L_p}\right)}$$

$$J_{top/bot}(L)=\frac{k_{top/bot}}{L_p^3}{\left(\frac{L_p}{L}\right)}^5{e}^{\left(2{\pi }^2\frac{L_p}{L}+\frac{L}{4L_p}\right)}$$

$$\text{centering}L=\frac{5}{\frac{1}{4L_p}+\frac{L_p{\theta }_a^2}{2}}\text{@add@centering}$$

Full text

I Introduction

Deformation of polymers is ubiquitous, elastic properties of these macromolecules are crucial for their dynamics. Biopolymers are abundant in nature and play vital roles in many biological processes Grosberg et al. (1997); Garcia et al. (2007); Bustamante et al. (2004); Nelson (1999), which not only depend upon the polymer structure, but also their physical properties Van de Velde and Kiekens (2002); Gosline et al. (2002); Kasas et al. (2004). Among biopolymers, DNA stands out as a case of its own. Understanding DNA deformation is crucial for the mechanistic grasp of vital cellular functions such as packaging of DNA compactly into viral capsids, the chromatin, formation of protein/DNA complexes and regulation of gene expression Garcia et al. (2007); Peters and Maher (2010). An all-important example of DNA deformation, relevant to a variety of biological processes that depend on its elastic properties, is DNA looping, which occurs in many prokaryotic Gralla (1991) and eukariotic Richmond and Davey (2003) systems. A number of regulatory proteins can loop DNA into various bent conformations, critical for regulation of many biological processes involving DNA Schleif (1992). Most notably, DNA is strongly bent in the nucleosome Kornberg (1974); Luger et al. (1997), which is the fundamental unit of genome packing: accessibility to genomic information in eukareotes is modulated by the strength of DNA-protein association Henikoff (2008); Fenley et al. (2018). Note that the majority of eukariotic genomic DNA (75-80%) is packed tightly into nucleosomes Richmond and Davey (2003). Nanostructures made directly of DNAMaffeo et al. (2016) or those that use DNA as a scaffoldMastroianni et al. (2009), can be influencedPark et al. (2008) by its mechanical properties on short length scales, providing yet another impetus to understand the strong bending regime of the DNA.

Experimental evidence on cyclization of DNA fragments shorter than $\sim$100 base-pairs points to the fact that strongly bent DNA – most relevant from a biological perspective – is considerably more flexible than expected from established models (worm-like chain) that work well within the weak bending regime. Yet, despite decades of experimental and theoretical effort, the story of how this arguably most important polymer behaves under deformation is far from complete, with controversies and new developments abound Vafabakhsh and Ha (2012); Vologodskii and Frank-Kamenetskii (2013); Savin et al. (2013); Zoli (2018); Sivak and Geissler (2012); Salari et al. (2015); Wu et al. (2015). The current state and the relevant terminology are briefly reviewed below.

The bending flexibility of a polymer is conventionally quantified in terms of its persistence length, $L_{p}$, a length scale below which the polymer behaves more or less like a rigid rod. Specifically, $L_{p}$ is defined as length of the polymer segment over which the time-averaged orientation of the polymer becomes uncorrelated; for fragments smaller than $L_{p}$, the thermal fluctuation alone is not enough to induce significant ($\sim$ 1 rad) bending Travers and Thompson (2004). Here we use this definition of $L_{p}$ to qualitatively separate the two bending regimes: if no significant bending is observed on length scales shorter than $L_{p}$, the polymer can be deemed weakly bent; otherwise the bending is assumed to be strong.

For double-stranded DNA, a variety of experimental techniques Baumann et al. (1997); Crothers et al. (1992, 1992); Vologodskaia and Vologodskii (2002); Shore et al. (1981); Hagerman (1988); Du et al. (2005), revealed that $L_{p}\approx 150$ bp or 500 Å. Based on the $L_{p}$ value and the above definition of strong bending, we conclude that most of the DNA in eukareotes is strongly bent. Indeed, since the nucleosome contains a stretch of double-stranded DNA of $\sim 150$ bp looped almost twice, the DNA in this complex can be considered as strongly bent.

Response of DNA to mechanical stress has been studied extensively  Frank-Kamenetskii et al. (1985); Shore et al. (1981); Marko and Siggia (1994); Baumann et al. (1997); Bloomfield (1997); Garcia et al. (2007); Lankaš et al. (2000); Cloutier and Widom (2004); Wiggins et al. (2006, 2005); Biswas et al. (2012); Bomble and Case (2008); Cloutier and Widom (2005); Du et al. (2005); Lankaš et al. (2006); Mazur (2007); Seol et al. (2007); Strauss and Maher (1994); Prévost et al. (2009); Peters and Maher (2010); Fields et al. (2013); Mathew-Fenn et al. (2008), leading to a consensus in modeling the weak bending regime. Arguably the most widely used simplified model of DNA bending is the worm-like chain (WLC) model. In the original WLC Kratky and Porod (1949); Du et al. (2005), the polymer is modeled as a continuous, isotropic elastic rod with its deformation energy being a quadratic function of the deformation angle. In the discrete version of WLC model, the bending energy of the polymer consisting of $N$ segments of length $l$ is given by:

[TABLE]

where $\theta_{i}$ is the angle between two consecutive segments (see inset of Fig. 1). While this simplistic model lacks some features of the real DNA, such as sequence dependence of its local mechanical properties, it nevertheless captures the key physics of weak polymer bending, which explains why the model is robust and is widely adopted to interpret experiment. Various theoretical models of DNA bending, including those that explicitly account for the sequence-dependence Coleman et al. (2003); Lankaš et al. (2003); Dixit et al. (2005); Fujii et al. (2007), were consequently developed that also assumed harmonic (quadratic) angular deformation energy of DNA. There is very little doubt that the Hookean, “elastic rod” models accurately describe many polymers in the weak bending regime Grosberg et al. (1997), including the double-stranded DNADu et al. (2005); Mazur (2007). Indeed, lowest order term of a Taylor series expansion of any well-behaved function around its local minimum is quadratic, which means that for small deviations from equilibrium, the response function can be considered harmonic. However, by the same logic it should be expected that beyond a certain threshold the bending energy may no longer be approximated by the quadratic term alone; investigations of possible influence of non-harmonic terms on the mechanical properties of double-stranded DNA is a relatively new area. Historically, only very large fragments (hundreds to thousands of base-pairs) were investigated Baumann et al. (1997); Crothers et al. (1992), which are well described by the WLC regardless of what happens on short length-scales Wiggins et al. (2006).

However, within the last decade or so, the prevailing view of DNA as a Hookean polymer was challenged by experiments that were able to investigate the flexibility of DNA on scales smaller than several $L_{p}$. Counter-intuitively, small DNA fragments ($\approx 100$bp) were found to have much higher probability of cyclization (spontaneous formation of loops) than that predicted by the WLC theory Cloutier and Widom (2005). This discovery sparked considerable controversy, which still remains unresolved. What is particularly puzzling is that strongly bent DNA appears less rigid than the DNA in the Hookean regime. Some of the follow-up experimental and theoretical work supports the validity of WLC even for tightly bent DNA Du et al. (2005); Vologodskii and Frank-Kamenetskii (2013); Mazur and Maaloum (2014), while others still show that short, tightly bent DNA is much more flexible Wiggins et al. (2006); Czapla et al. (2006); Vafabakhsh and Ha (2012) than previously thought, in a manner that can not be described by a harmonic model Wiggins et al. (2006).

Several theoretical models have been proposed to account for the unexpectedly high flexibility of strongly bent double-stranded DNA. One popular modelYan and Marko (2004) – the meltable WLC or MWLC – postulated that the extra flexibility stems from formation of small local “bubbles” of single stranded DNA, which is much softer than the double helix. However, the degree of softening provided by the mechanism alone was later foundForties et al. (2009) to be inadequate to fully explain the very sharp bends in DNA observed experimentally; in atomistic simulations, negative super-coiling was required to induce such bubbles in DNA mini-circlesMitchell et al. (2011). An early modelCrick and Klug (1975), put forward well before the unusual DNA flexibility was discovered experimentally, suggested that the energy of a bent double-helix could be lowered by formation of sharp, $\sim 90^{o}$ kinks that maintain the Watson-Crick pairing along the helix. Sharp kinks were indeed observed in a pioneering atomistic simulationLankaš et al. (2006) some thirty years later, but subsequent improvements in the simulation methodology indicated that these were only induced at a high bend angle equivalent to those occurring in circles of just $45$ base-pairsCuruksu et al. (2009), while experimental softening of the DNA is seen experimentally for circles as large as $\sim 106$ base-pairsVafabakhsh and Ha (2012). Sharp kinks in double-stranded DNA can be introduced empirically into the WLC model, e.g. by adding freely-bending hinge elements to the WLC chain, leading to a kinkable WLC, or KWLCWiggins et al. (2005). A non-linear empirical bending potential that allows for the possibility of $\sim 90^{o}$ kinks in double-stranded DNA was recently proposedVologodskii and Frank-Kamenetskii (2013), but its physical origins, the critical value of the DNA curvature at which the kink occurs, and the corresponding energy gain remained unknownVologodskii and Frank-Kamenetskii (2013). At the same time, a purely linear empirical bending potential was shownWiggins et al. (2006) to describe the softer DNA seen in AFM experiments, although the origin of the linear regime and its parameters (e.g. critical bend angle where the linear regime begins ) remained unclear. Are the kinking and the linear regime just two manifestations of a deeper underlying principle?

In summary, the nature of the effective bending energy of double-stranded DNA in the strong bending regime as well as the precise connection to the observed softening of the polymer is not fully clear. The influence of mechanical constraints on this connection remains unexplored. It is unclear how the softening of strongly bent DNA stems from its atomic-level structure and interactions. From a more philosophical standpoint, is hard to believe that very special models are needed to describe the bending of the DNA; rather it is more like that the curious case of the DNA is just a special case of a broader underlying theory applicable to all polymers.

In this work we propose, and verify against available experiment, a unified theoretical description of polymer bending that treats the weak and strong bending regimes on the same footing, guided by a simple physical principle. The proposed framework does not rely on ad-hoc postulates; instead, it shows how the apparent softening of strongly bent DNA follows naturally from a specific mathematical property of the experimentally-derived bending energy. Simulations suggest an atomistic explanation for the specific shape of the bending energy function.

II Methods

II.1 DNA bending energy from experimental data

A statistically significant, diverse set of several hundred PDB structures of protein-DNA complexes was investigated previously in Ref. Du et al. (2005). The probability distribution of the experimental DNA bending angles was used in Ref. Du et al. (2005) to approximate the bending energy $E(\theta)$ (per base pair) as a fourth oder polynomial: $E(\theta)=203.1\theta^{2}-552.7\theta^{3}+416.8\theta^{4}$ (where $\theta$ is in radians and $E(\theta)$ is in units of $kT$). Here we use this $E(\theta)$ to represent the experimental effective bending energy of the double-stranded DNA, blue line in Fig. 2.

II.2 Atomistic MD simulations of closed DNA loops

To avoid “end effects”, and make a close connection with DNA cyclization experiments, we employed closed DNA circles to estimate their effective bending energy $E(\theta)$ per base pair. DNA circles of various sizes (50-400 bp) were generated using NABMacke and Case (1998) (AmberTools) for sequence poly(dA).poly(dT), helical repeat of 10 bp, and other parameters of B-DNA as specified in NAB. We deliberately chose this simple, uniform sequence to focus on the basic physics of DNA deformation.

All of the atomistic MD simulations were performed within AMBER-12 package, using ff99bsc0 force-field. The Generalized Born (GB-HCT, AMBER option igb=1) implicit solvation model was used to treat solvation effects, including 0.145M of monovalent salt. No long-range cut-off was employed. The model’s performance in atomistic simulations of DNA, including studies of its deformationBomble and Case (2008), is well establishedOnufriev (2010). Two critical advantages of the implicit solvation over the more traditional explicit solvationOnufriev and Izadi (2018) made the former the method of choice in this work. These are the superior simulation efficiency for large DNA structuresAnandakrishnan et al. (2015) and the straightforward manner in which their energies, including free energy of solvent re-arrangement, can be estimatedOnufriev (2010) within the implicit solvation framework.

All DNA circles were initially minimized for 1000 steps with “P” atoms restrained their original positions with a force constant of 1.0 kcal/mol/Å2 to enforce the circular shape. Each system was then heated to 300K and equilibrated for 100 ps with the same restraints as for the minimization. Shake was used to constrain the hydrogen atoms; we employed 2 fs time-step for the atomistic simulations. Finally, we generated 1 ns long MD trajectory for each circle at 300K, with “P” atoms also restrained with a force constant of 0.1 kcal/mol/Å2, sufficient to support the near perfect circular shape of the fragment, but allowing for local re-arrangements. The energies and their components, including the electrostatic, VDW, bond, etc. were saved every 20 fs, and averaged over the whole trajectory. The relatively short simulation time allowed us to simulate even the largest of the circles; it is justified by the use of the strong positional restraints, which permit only local, very fast structural re-arrangements. For smaller circles we verified that increasing the simulation time by an order of magnitude had negligible effect on the computed averages. In Fig. 4, the energy per bp was computed as the difference between per bp potential energies of the given circle and the largest circle simulated, which is virtually unbent.

II.3 Coarse-grained simulations of DNA loops

ESPResSo Limbach et al. (2006) was used to create and simulate coarse-grained closed loops of DNA of different sizes, from 6 to 600 bp long. A single bead of the appropriate mass represents one base-pair of B-DNA; the bead-bead distance was set to $3.3$ Å , corresponding to the average distance between base pairs in canonical DNA. The bonds between the beads were made virtually inextensible (very large coefficient of the quadratic bond stretching energy); the bond angle potential (effective bending energy) between neighboring beads was defined to have the same form as in Fig. 2, that is correspond to the bending potential inferred from the experimental dataDu et al. (2005). No further bead-bead interactions or constraints on the loop geometry were imposed. The loops were simulated at $T=300K$, and energy-minimized using steepest descent.

II.4 Coarse-grained simulations of confined DNA fragments

Protein-DNA complex.

ESPResSo Limbach et al. (2006) was used to create and simulate a 20 bead long fragment of “DNA” bound to a spherical charged ”protein”, Fig. 5. The beads and their interactions were set up as described above, with the following modifications. The end beads were not linked to create a loop. Each bead carried a unit charge $q_{s}=$-1; The bead charges interacted only with a positive charge $Q$ of the “protein”, represented by a spherical impenetrable constraint of radius $R$. In addition, two impenetrable walls were placed above and below the charge $Q$ to minimize out-of-plane bending of the “DNA”. The confining charge $Q$ was varied from 10 to 1000, effectively sampling two orders of magnitude of confinement strength (defined here as $|Q/q_{s}|$). The constraint radius $R$ was also varied to sample various curvature values of the “protein”, and thus various total bending angles of the confined “DNA”, Fig. 5.

A nucleosome model.

For the nucleosome model, the system described above was modified to mimic the confinement of DNA around realistic histone core. The DNA fragment size was increased to 147 bp, and the non-bonded interactions between monomers were turned on for an additional realismFenley et al. (2010). The fragment was confined around a cylinder of fixed diameter $R=\sim$100Å, and the walls were placed $\sim$50Å apart (approximate dimensions of the nucleosome complexFenley et al. (2010)).

III Results and Discussion

III.1 The proposed unified framework of polymer bending

We begin with a useful analogy from classical thermodynamics that connects system’s stability to convex properties of its governing potential. For example, for a system to be stable against a macroscopic fluctuation in energy, the entropy of the system as a function of energy, $S(E)$, must be concave (non-convex). Any chord connecting two points on a graph of $S(E)$ must lie below the curve itself in order to satisfy the second law of thermodynamics (maximum $S$). Conversely, the inverse function $E(S)$ must be convex. If, however, $E(S)$ is not convex over some region, the system phase-separates once this region is reached, with the properties of the two phases corresponding to the end points of the convex hull of the non-convex region. The actual, physical average energy of the system follows the convex hull, which makes the energy manifestly convex. This very general reasoning, with appropriate choice of the perturbation coordinate and potential, is applicable to phase transition of single species polymers (Flory-Huggins TheoryFlory (1953)), as well as to stretching of polymersSavin et al. (2013) and other materialsBabicheva et al. (2013). Here we use the analogy to develop a general framework that describes polymer response to bending, weak and strong, on the same footing.

Consider a polymer chain made of $N\gg 1$ inextensible, identical monomer segments with effective bending deformation energy $E(\theta_{i})$ for each bending site, where $\theta_{i}$ is the angle between two successive segments (see inset of Fig. 1). Here we assume that the effective $E(\theta)$ takes into account all the interactions, short- and long- range, between the monomers. For notational simplicity, in what follows we ignore the difference between $N$ and $N-1$ for large $N$. The total energy of the polymer is $E_{chain}=\sum_{i}^{N}E(\theta_{i})$, and without loss of generality we assume no intrinsic bends, i.e. $E(0)=0$. Just like in WLC, we assume isotropic bending energy, which is a reasonable assumption for DNA fragments longer than 2 helical repeats or 20 bpVologodskii and Frank-Kamenetskii (2013). For the moment, we further assume no torsional degrees of freedom. In order to induce an average non-zero bend in the chain, the polymer must be constrained, and the problem of finding the equilibrium polymer conformation is reduced to minimizing $E_{chain}$, subject to the specific constraint of the problem. Here we assume that entropic effects are relatively small at length scales of interest ($\lesssim L_{p}$) – an assumption that we explicitly confirm below by numerical experiments.

We begin by considering a very special case of a uniformly bent polymer – constrained to have the same constant curvature along the entire chain. By construction, such a polymer consists of identically bent segments with each bending angle $\theta_{i}$ equal to the average deformation angle, $\bar{\theta}=N^{-1}\sum_{i}^{N}\theta_{i}$, and its total energy is $NE(\bar{\theta})$.

Next, consider a more realistic situation where the polymer bending is enforced by a much less restrictive constraint: that the sum of the bend angles between the monomers remains constant, $\alpha=\sum_{i}^{N}\theta_{i}=const$. Note that this constraint alone does not fully define the geometry of the polymer. A closed planar loop, with the first and last segments linked, is a relevant example for which the constraint is satisfied; $\sum_{i}^{N}\theta_{i}=2\pi$, from elementary geometry of polygons, see also the SI. Mathematically, the problem of finding the minimum energy conformation of the polymer is that of energy minimization under the specific constraint:

[TABLE]

where we make a clear distinction between $\mathbfcal E$, which is the average bending energy per bending site in the minimum energy state of the polymer, and $E$ corresponding to the uniform bending. Using Lagrange multipliers, Eq. 2 can be reduced to $\min\{E(\theta_{1})+\dots+E(\theta_{N})-\lambda(\theta_{1}+\dots+\theta_{N}-\alpha)\}$. Differentiating with respect to $\theta_{i}$ gives a set of equations $\partial_{\theta_{i}}E(\theta_{i})-\lambda=0$ (for all $i$) which leads to a set of equalities $\partial_{\theta_{1}}E(\theta_{1})=\partial_{\theta_{2}}E(\theta_{2})=\dots=\partial_{\theta_{N}}E(\theta_{N})$. For a convex functional form of $E(\theta)$, $\partial_{\theta}E(\theta)$ monotonically increases with $\theta$, and therefore the equalities are satisfied only if $\theta_{1}=\theta_{2}=\dots=\theta_{N}$: the polymer is always uniformly bent, that is each segment is bent through the same angle $\theta=\bar{\theta}$ and $\mathbfcal E(\bar{\theta})=E(\bar{\theta})$. However, for a non-convex function such as one shown in Fig. 1, there can be more than one value of $\theta$ that satisfies $\partial_{\theta_{1}}E(\theta_{1})=\partial_{\theta_{2}}E(\theta_{2})=\dots=\partial_{\theta_{N}}E(\theta_{N})$: $\partial_{\theta_{i}}E(\theta_{a})=\partial_{\theta_{i}}E(\theta_{b})$ for some $\theta_{a}<\theta_{b}$. Of special importance are $\theta_{a}$ and $\theta_{b}$ that mark the beginning and the end of the convex hull of $E(\theta)$ – the segment of a straight line tangent to the non-convex function at two points, such that for any argument between these two points the value of the function at the argument is greater than that of the convex hull line at the same argument, Fig. 1. One can demonstrate, see SI, that for bend angles $\bar{\theta}$ in the convex hull interval, $\theta_{a}<\bar{\theta}<\theta_{b}$, a uniformly bent chain is no longer the stable minimum energy conformation of the polymer. Instead, the stable minimum is achieved when the distribution of bend angles is bi-modal: each segment is bent through one of the two bending angles $\theta_{a}$ or $\theta_{b}$. This general point is illustrated in SI for a model polymer chain described by a non-convex bending potential relevant to the case of the DNA.

In what follows we derive an explicit expression for $\mathbfcal E(\theta)$ for $\theta_{a}<\bar{\theta}<\theta_{b}$. In the minimum energy conformation, let $0<p<1$ represent the fraction of all the bending sites that are in the state $\theta_{b}$ and $1-p$ the fraction of the remaining sites in the state $\theta_{a}$. The total bending angle in terms of $\theta_{a}$ and $\theta_{b}$ is then given by $Np\theta_{b}+N(1-p)\theta_{a}=N\bar{\theta}=\alpha$, and the bending energy per monomer in the non-convex region is $\mathbfcal E(\bar{\theta})=pE(\theta_{b})+(1-p)E(\theta_{a})$. Rewriting $p=(\bar{\theta}-\theta_{a})/(\theta_{b}-\theta_{a})$, we arrive at

[TABLE]

Therefore, in the non-convex region, the actual polymer energy per bending site, $\mathbfcal E(\bar{\theta})$ corresponding to the stable minimum energy state, is a linear function of the average deformation $\bar{\theta}$. Clearly, $\mathbfcal E(\bar{\theta})<E(\bar{\theta})$ within the convex hull interval, Fig. 1.

To arrive at a general theory that can account for both the weak and strong bending regimes simultaneously, we use the form of Eq. 3 for the strong bending regime, while retaining WLC for the weak bending. In the proposed Energy Convex Hull(ECH ) model, the average per segment (e.g. per base-pair) bending energy is described by an everywhere differentiable piece-wise polynomial function: quadratic WLC (Eq. 1) for $\bar{\theta}<\theta_{a}$, and a linear function – convex hull of $E(\theta)$ – for $\theta_{a}<\bar{\theta}<\theta_{b}$:

[TABLE]

where $L_{p}$ is the accepted persistence length, well established for the weak bending regime; here it dimensionless, expressed in terms of the number of bending sites (e.g. number of base pairs for DNA loops). In this work we are not interested in the extreme strong bending regime $\bar{\theta}>\theta_{b}$, since for the DNA this regime would correspond to loops smaller than 10 bp. Such small loops are likely physically impossible due to steric constraints, and are much smaller than those observed in cyclization studiesVafabakhsh and Ha (2012); Lionberger et al. (2011). Thus, the only key parameter that ECH theory inherits from the input effective bending energy, $E(\theta)$ in Fig. 2, is the value of $\theta_{a}$, which enhances robustness of the theory to inevitable imperfectionsDu et al. (2005) of the input bending energy profile. For example, a uniform re-scaling $E(\theta)\rightarrow\lambda E(\theta)$ would leave the x-coordinates $\theta_{a}$ and $\theta_{b}$ of the convex hull double-tangent segment unchanged because the derivatives would be re-scaled by the same $\lambda$. Further discussion of the robustness of ECH model to its parameters can be found below and in SI.

III.2 Bending of a circular loop, weak and strong

While many different types of constraints can be physically realized, one of the most important ones is the closed loop constraint, which is also used in DNA cyclization experimentsShore et al. (1981); Du et al. (2005); Cloutier and Widom (2005) criticalWiggins et al. (2005) for investigating the strong bending regime. Consider the case of a single closed loop $\alpha=\sum_{i}^{N}\theta_{i}=2\pi$. From Eq. 4, the total bending energy of a closed loop of total length $L$ (number of base pairs, corresponding to “N” in Eq. 2 ) is given by $\mathbfcal E_{loop}=\mathbfcal E(\bar{\theta})L$. Since $\bar{\theta}=\frac{\alpha}{L}=\frac{2\pi}{L}$, the bending energy of the loop is:

[TABLE]

Note that, where defined, the new function $\mathbfcal E_{loop}$ depends on just one new parameter: $\theta_{a}$ – lower boundary of the non-convex domain. Although we tacitly assumed the loop to be confined to a 2D plane to simplify the derivations, our unconstrained coarse-grained simulations of closed loops at 300K demonstrate, Fig. 2, that the assumption has little effect on our key conclusions.

III.3 Application to double-stranded DNA

The preceding discussion was not restricted to the case of DNA: non-uniform, two-phase bending, and the corresponding linear bending regime can be a feature of any polymer. However, since DNA is arguably the most important polymer, and it exhibits looping in many different biological systems, we will focus on double-stranded DNA for the rest of the study. An effective bending energy (per bp) calculated from a statistical analysis of experimental PDB structures of DNA-protein complexes Du et al. (2005) is shown in Fig. 2. This effective bending energy function has a non-convex region, and thus a convex hull, the end points of which are $\theta_{a}=2.2^{\circ}$ and $\theta_{b}=35.8^{\circ}$, corresponding to fragment lengths of $L\sim 160$ and $\sim 10$ bp respectively for DNA closed loops.

Coarse-grain molecular dynamics simulations at 300K (see Supplementary Material) demonstrate that polymers with this effective bending energy between monomers exhibit all of the key features discussed above. For large loop sizes the bending angles are small (weak bending) – the system samples the convex (harmonic) region of the energy function, Fig. 2, and the distribution of bend angles is uni-modal. However, as the loop size decreases, the average angle per bending site $\bar{\theta}$ increases, eventually crossing the $\theta_{a}$ threshold. Once this happens, the energy of the system per bending site increases linearly with $\bar{\theta}$, and the distribution of bend angles becomes bi-modal, until the system reaches the upper boundary of the convex hull at $\theta_{b}$.

III.4 Comparison with DNA cyclization experiments

Most experimental cyclization results are expressed Cloutier and Widom (2005); Du et al. (2005) in terms of the Jacobson-Stockmayer j-factor, which estimates the probability that a linear polymer of length $L$ forms a closed loop by joining its cohesive ends Shore et al. (1981); Jacobson and Stockmayer (1950). While Monte-Carlo based numerical approaches to compute j-factor existSivak and Geissler (2012), here we use a well-established Shimada and Yamakawa (1984); Allemand et al. (2006) analytical expression for the j-factor of an unconstrained closed loop:

[TABLE]

where $\frac{1}{L_{p}^{3}}\left(\frac{L_{p}}{L}\right)^{5}\exp\left(\frac{L}{4L_{p}}\right)$ accounts for the entropic contribution, averaged over possible looping geometries, and $\exp\left(-\frac{\mathbfcal E_{loop}}{k_{B}T}\right)$ is the energy penalty of bending the DNA fragment to form the loop. We note that the $k$ depends, in a complex manner, on the loop closing geometry and can be expected to remain invariant over a relatively short range of loop lengths $L$, within the same experiment.

To make a direct connection with cyclization experiments for non-integer numbers of helical repeats, we modulate the torsionally independent loop energy from Eq. 6 with $cos(2\pi L/h)$, where we assumed the helical repeat $h=10$ bp per turn. The agreement with the cyclization experiment is robust with respect to the precise value of the helical repeat, see SI. This form of the modulating factor is adopted from Ref. Shimada and Yamakawa (1984) to account for the periodic variation of the j-factor due to the torsional component of the energy Shimada and Yamakawa (1984). This simple way of accounting for non-integer numbers of helical repeats is sufficient for the purpose of testing key predictions of ECH vs. WLC, and does not affect the comparison with the over-all (envelope, average) behaviorWiggins et al. (2005) of experimental j-factors, see also Table S1 in SI. We use $\mathbfcal E_{loop}(L)$ defined in Eq. 5 for ECH and $\mathbfcal E_{loop}(L)=2\pi^{2}k_{B}T\frac{L_{p}}{L}$ for all $L$ in the case of WLC. The proposed ECH model and WLC are compared with the most recent experiment Vafabakhsh and Ha (2012) in Fig. 3.

As seen from Fig. 3, ECH  leads to an excellent agreement with the cyclization experiment, while the $j$-factors predicted by conventional WLC are off by several orders of magnitude in the strong bending regime (WLC is known to work well in the weak bending regime where it coincides with ECH by construction). The agreement of ECH with the experiment is robust to the value of its key input parameter $\theta_{a}$, see below and SI.

Cyclization of very short loops.

Counter-intuitively, the predicted envelope function for ECH j-factor, which is essentially Eq. 6, has a minimum near 45 bp and begins to increase for even smaller loops, whereas for WLC j-factor decreases sharply for small loops. This completely counter-intuitive behavior of the cyclization probability for very tight loops predicted by ECH is borne out by experiment, Fig. 3; its physical origin is explained below. The two experimental points at L=50 and L=40, which support the counter-intuitive prediction of the theory, were not available to us until after the ECH framework was fully developed and tested against published dataVafabakhsh and Ha (2012) for larger circles.

The over-all variation of the j-factor as a function of the loop length for both models is governed by the interplay between the entropic and the mechanical bending energy costs $\mathbfcal E_{loop}(L)$ of forming the loop. For small loops, the entropic penalty of forming the loop decreases with the loop size $L$; however, $\mathbfcal E_{loop}(L)\rightarrow\infty$ for small $L$ within WLC, which leads to a steep decrease in the over-all cyclization probability. In contrast, ECH loop energy, Eq. 5, approaches a constant for $L\rightarrow 0$, which explains why the corresponding j-factor reaches a minimum and then begins to increase for small enough $L$, Fig. 3. This very different qualitative behavior of WLC and ECH j-factors for small loops can be used as a discriminating experimental test of the models. The predicted minimum value of the j-factor can be used to further discriminate between models that exhibit the minimum: for example, both KWLCWiggins et al. (2005) and a recent versionSivak and Geissler (2012) of MWLC predict the minimum, but the loop sizes at which the minima occur are substantially different from the $\sim 45$ base pairs predicted by ECH .

Within the proposed ECH framework of polymer bending, the central role is played by convexity properties of the effective bending energy between individual monomers. For the DNA, we used the energy profile inferred from statistical analysis of experimental structures of protein-DNA complexes (Fig. 2) – the energy has a clear non-convex region, responsible for the “softer”, linear bending mode of short DNA loops. The same general considerations will hold for any effective bending energy that has a distinct non-convex region regardless of its originSalari et al. (2015), including a kinkable WLC (kWLC) potentialVologodskii and Frank-Kamenetskii (2013). Thus, even though ECH explains experimental results perceived to be in contradiction with WLC, there is no fundamental contradiction between the new framework and the conceptual basis of WLC.

III.5 Origin of the non-convex bending energy of DNA.

To investigate, qualitatively, the physical origin of the non-convexity of the DNA effective bending energy we employed all-atom Molecular Dynamic (MD) simulations of uniformly bent DNA circles of a wide range of sizes, from small to very large, corresponding to almost unbent DNA, see “Methods”. Specifically, we examined the average bending energy per base pare. The total bending energy profile obtained from these simulations, along with the breakdown into components of different physical origin, are shown in Fig. 4; one can clearly see a prominent non-convex region, in qualitative agreement with the experiment, Fig. 2. The key parameter $\theta_{a}\approx 1.5^{\circ}$ from the MD simulations, which is not all that different from the value of $2.2^{\circ}$ inferred from the experimental data, Fig. 2. Some discrepancy is likely due to sequence effectsLavery et al. (2010); Wang et al. (2017), force-field issuesCheatham and Case (2013); Savelyev and MacKerell Jr (2014), or the fact that the experiment-based potential in Fig. 2 may itself deviate from reality to some extent, as noted in the original publicationDu et al. (2005). Importantly, the use of MD-derived $\theta_{a}=1.5^{\circ}$ in Eqs. 5 and 6 results, see SI, in virtually the same close agreement with the cyclization experiment we have seen Fig. 3, which is based on $\theta_{a}=2.2^{\circ}$ derived from experiment. This insensitivity of the prediction of ECH to the value of its key input parameter points again to the robustness of the framework. The following qualitative conclusions can be made from the MD-based analysis of the DNA bending, Fig. 2. For small bending angles, the total energy is reasonably well approximated by a quadratic function. However, once the bending reaches the transition angle $\theta_{a}$, the VDW energy decreases at a rate faster than the increase of the other terms combined, which results in a non-convex region of the total $E(\theta)$, Fig. 4. It is this sharp decrease in the VDW contribution that gives rise to the existence of a non-convex region in the DNA bending energy. Further analysis, (inset in Fig. 4), reveals that it is the attractive component of the VDW interactions between DNA backbone atoms (backbone-backbone), rather than base stacking, that is critical to the counter-intuitive sharp decrease in the total bending energy, see Supplementary Material for further atomistic details. The key role of the backbone-backbone VDW term suggests that it is the overall structure of DNA, rather than sequence details, that is responsible for the existence of the convex hull in the polymer’s effective bending profile; variations in the DNA sequence may alter the range of bending angles over which the convex hull exists.

It is worth mentioning that local “bubbles” of broken WC bonds do not occur in our atomistic MD simulations of DNA circles where the uniform bending is enforced by constraints on the phosphorous atoms, see “Methods”. Yet, these simulations yield a non-convex profile of the DNA bending energy, Fig. 4), which, as we have demonstrated, always leads to the existence of linear “soft” bending regime. Thus, the simulations suggest that local DNA melting may not be necessary to explain the high flexibility of strongly bent DNA and the stark deviation of experimental j-factors from WLC predictions. We stress that we do not rule out “bubbles” of broken WC bonds in actual sharply bent DNA; instead, we predict that if WC bond breaking were suppressed experimentally, the *qualitative * picture of sharply bent DNA being much softer would remain, with the experimental j-factors still deviating from the WLC in a way qualitatively similar111Quantitative details can be different if WC bond breaking is suppressed. Note that the effective loop bending energy of ECH theory in Fig. 3 comes from a statistical analysis of real protein-DNA complexes. Consequently, the ECH effective energy with parameters used in that figure, $\theta_{a}=2.2^{\circ}$ and $\theta_{b}=35.8^{\circ}$, implicitly accounts for broken WC bonds if these occur in the DNA of the complexes. to what is currently observed in experiment, Fig. 3.

An analogy can be made here with the physics behind the DNA overstretching plateauSmith et al. (1996); van Mameren et al. (2009), where the polymer extension occurs at constant force, and the stretching energy grows linearly with the polymer extension. This peculiar regime can be explainedSavin et al. (2013) via the same main argument used in the current work – the existence of a non-convex region in the polymer deformation energy. In the case of DNA overstretching, experiments have demonstrated convincinglyPaik and Perkins (2011) that WC bond breaking is not required for the existence of the characteristic plateau on the force-extension diagram.

III.6 Beyond closed loops: a protein-DNA “complex”.

The proposed framework is based on one main assumption: despite constraints, the polymer chain is still free to explore sufficient conformational space to search for minimum energy. So far, we focused on DNA loops because of direct connection to key cyclization experiments; the single constraint $\sum_{i}^{N}\theta_{i}=\alpha=2\pi$ is minimally restrictive. However, other realistic scenarios of DNA bending, notably in protein-DNA complexes, may involve different types of constraints that can confine the polymer strongly enough to potentially violate the main assumption to various degrees. Here we investigate to which extent our main conclusion – deformation energy of strongly bent DNA follows the convex hull of $E(\theta)$ – may still hold in a model of protein-DNA complex, Fig. 5 and “Methods”. We vary the total positive charge $Q$ of the cylindrical core to modulate the electrostatic attraction of the negatively charged polymer (monomer charge $q_{s}$) to the core surface of the “protein”, and, hence, the degree of the polymer confinement. In the limit of very strong confinement ($|Q/q_{s}|\rightarrow\infty$), the polymer is forced to be confined to a circular, uniformly bent path on the surface of the cylindrical core, and has very few degrees of freedom left to explore in this regime, solid red line in Fig. 5. The average bending energy in this case follows the given functional form of $E(\theta)$ (red dashed line in Fig. 5), and ECH clearly does not apply. As we decrease the confinement strength, the polymer is allowed to assume non-uniform bending conformations while lowering its total bending energy. The effective bending energy per monomer begins to approach the convex hull (solid green and blue lines), making ECH more applicable. In the case of the weakest confinement (solid purple line), the polymer is still loosely bound to the core, but is allowed to relax almost completely. This is the limiting case described by our ECH model: the resulting energy per monomer follows the convex hull fairly closely

We argue that it is this low confinement regime, where ECH is relevant, that describes the real nucleosome Luger et al. (1997) – arguably the most important DNA-protein complex. To illustrate, we model a “variable confinement nucleosome” by a coarse-grained 147-bp DNA fragment placed next to a cylinder with relative dimensions of the actual histone core Fenley et al. (2010), see Methods; as the core charge $Q$ is increased, the whole fragment starts to wrap around the cylinder once the confinement strength is $|Q/q_{s}|\geq 90$. At this value of the DNA confinement, the energy cost of pulling away a fragment of $\sim 20$ bp in our model is $\approx 10k_{B}T$, which is comparable to $\approx 6k_{B}T$ estimated from experiment Garcia et al. (2007) for the fragment of the same length in the actual nucleosome. Moreover, even for higher degrees of confinement, up to $|Q/q_{s}|\simeq 200$ ($\simeq 20k_{B}T$ to pull away a 20 bp fragment), the corresponding blue lines in Fig. 5 still approximate the convex hull, and so ECH is still likely applicable, at least qualitatively.

IV Conclusion

It is now well established that slightly bent DNA behaves like an elastic rod – the deformation energy is a quadratic (harmonic) function of the deformation. However, recent experiments demonstrated that strong bending of small DNA fragments could no longer be described within this classical model.

Here we have proposed a novel framework for bending of polymers, which is based on the consideration of convex properties of the effective bending energy between successive monomers. Within the framework, the bending energy is harmonic for small bends, but once the average deformation reaches the convex hull of the effective bending energy function, a “phase transition” to the strongly bent regime occurs, in which the system’s energy is a linear function of the average bending angle. In this regime, which persists for as long as the average deformation is within the convex hull interval, the two states of bending co-exits: some segments are bent weakly, while others are bent strongly, with the proportion of the latter increasing with the increased average bend (e.g. shorter loops). The transition point from the harmonic to the linear bending regime occurs at the beginning of the convex hull segment – this point plays a special role in the new theory. These general considerations are expected to hold for any polymer with an effective bending energy that has a distinct non-convex region, regardless of its originSalari et al. (2015); Vologodskii and Frank-Kamenetskii (2013).

For generic “sequence-averaged” double-stranded DNA considered here, we conclude that the effective bending deformation energy becomes non-convex for strong bends greater than $\sim 2^{\circ}$/bp, which corresponds to circular loops shorter than $\sim 160$ bp. The conclusion about the DNA bending energy being non-convex relies on an analysis of a large number of experimental protein-DNA complexes, and is consistent with the shape of the bending energy inferred from atomistic MD simulations. The simulations also yield a qualitatively similar value for the bend angle that marks the onset of the linear bending regime. Further, atomistic simulations of DNA circles reveal that the attractive short-range Lennard-Jones interactions between the backbone atoms are key for the underlying non-convexity of the DNA effective bending energy, leading to the linear bending regime. We use MD simulations only for general reasoning, which is robust to details of the simulation protocol.

In this work our focus is the main principle; future refinements of the ECH theory may be able to account for details not considered here, such as sequence dependence of the DNA bending energy, the influence of torsional stress and supercoiling, etc. We have also just barely touched upon structural consequences of ECH , such as the number and distribution of “kinks” in tightly bent DNA. An analysis of these features will likely lead to more verifiable predictions of the theory. Likewise, we have derived specific mathematical expressions for bending under only one type of constraint; other relevant types of constraints need to be considered in more detail to complete the theory. Based on our analysis, the key conceptual features of ECH will likely hold.

The new theory does not contradict the conceptual basis of the classical models of DNA bending such as WLC, but also agrees with recent experimental cyclization data on strongly bent small DNA circlesVafabakhsh and Ha (2012). A completely counter-intuitive prediction that cyclization probability reaches a minimum for very small loops has proved to be consistent with additional experimental data points, not available to us when we made the prediction.

We believe that the novel general framework can be used to analyze, at least conceptually, many other scenarios of strong polymer bending, and should help interpret future experimental observations.

V Funding

This work was supported in part by the National Institutes of Health [R01 GM099450] and the National Science Foundation [MCB-1715207].

VI Acknowledgments

We thank Taekjip Ha for sharing with us unpublished experimental data. We thank Igor Tolokh for many detailed and insightful comments.

The authors acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper.

VII Author Contributions

AD, AM and AVO performed the research and wrote the manuscript. AVO designed the research.

**SUPPLEMENTARY MATERIAL **

VIII Convex vs. non-convex bending energy functions

Below is an argument for the difference in bending behavior of polymers characterized by a convex vs. non-convex bending energy. As in the main text, we assume that the polymer deformation energy has only the bending component, which is isotropic. In all cases, bending is induced by constraint $\sum_{i}^{N}\theta_{i}=N\bar{\theta}=const$, where $\bar{\theta}$ is the average bending angle of the polymer chain made of $N$ segments.

Without the loss of generality, consider the uniformly bent conformation of a polymer with $N=2$, and let’s investigate its stability to perturbation. A perturbation $\Delta\theta$ that reduces the bend angle at one site means that the other bending site must increase its bend angle by $\Delta\theta$ in order to satisfy the constraint $N\bar{\theta}=const$. This perturbation changes the total energy $E_{chain}$ by $E(\bar{\theta}+\Delta\theta)+E(\bar{\theta}-\Delta\theta)-2E(\bar{\theta})$. If $E(\theta)$ is a convex function (e.g. the black curve in Fig. 6), the perturbed system will have a higher energy than in the initial uniformly bent case. This is because, by definition, a convex curve always lies below its chords, so that $2E(\bar{\theta})<E(\bar{\theta}+\Delta\theta)+E(\bar{\theta}-\Delta\theta)$. Thus, the perturbation leads to increase in system energy which implies that the system was at a stable equilibrium. Therefore, the minimum energy conformation of a polymer with a convex effective bending energy is always that of a uniformly bent chain.

However, if the function $E(\theta)$ has a non-convex region, e.g. the blue curve in Fig. 6, then for any $\bar{\theta}$ in that region, and $\Delta\theta$ that does not take the system outside of it, $2E(\bar{\theta})>E(\bar{\theta}+\Delta\theta)+E(\bar{\theta}-\Delta\theta)$ , which means it is possible to lower the energy of the polymer further by non-uniform bending. Namely, one site is now bent through $\bar{\theta}-\Delta\theta$, and the other through $\bar{\theta}+\Delta\theta$, with the new value of the average chain energy per site, $\frac{1}{2}E_{chain}$ falling on the midpoint of the line (dashed red line in Fig. 6) connecting the two new bending states on the energy curve. The longer the cord connecting the two perturbed states at $\bar{\theta}-\Delta\theta$ and $\bar{\theta}+\Delta\theta$, the larger the energy gain $2E(\bar{\theta})-E(\bar{\theta}+\Delta\theta)-E(\bar{\theta}-\Delta\theta)$ due to the non-uniform bending, for as long as the cord is completely below the $E(\theta)$ curve. The largest energy gain, and thus lowest possible $E_{chain}$, is achieved for the limiting cord that is the convex hull of $E(\theta)$ – the line segment tangent to the non-convex function at two points, such that between these two points the value of the function is greater than that at any point of the line segment. For this limiting case, $\bar{\theta}-\Delta\theta=\theta_{a}$, $\bar{\theta}+\Delta\theta=\theta_{b}$.

IX Concave vs. convex closed loops.

The exterior angle sum theorem is valid for convex polygons, directly applicable to the sum of the tangent vectors of a convex closed curve. We do not consider non-convex (concave) closed curves here because these can not minimize the deformation energy of a closed inextensible elastic loop, at least as long as the effective bending energy $E(\theta)$ is a monotonic function of $\theta$, or curvature $\kappa$. (The existence of a non-convex region in $E(\theta)$ does not imply a non-monotonic $E(\theta)$. )

The outline of an intuitive proof idea is as follows. Any concave curve can be transformed – “banged out” – to eliminate a local non-convex region while bringing the bending energy down in the process, see Fig. 7 for an illustration of the process for a shallow “dimple”. For a deeper “dimple”, the procedure may involve two steps. First, the entire convex portion of the curve is “stretched out” via a uniform re-scaling of its curvature $\kappa(s)\rightarrow\lambda\kappa(s),0<\lambda<1$. The bending energy of the convex portion of the curve will become lower during this step. Since the curve is inextensible, and there are no breaks, points A and B will move further apart as a result, thus also reducing the curvature of the non-convex “dimple” that spans the $|AB|$ segment, and hence lowering its bending energy. After that, the step in Fig. 7 is applied, making the resulting curve convex, and lowering its bending energy further. We do not pursue a more formal proof based on variational calculus.

X Non-convex bending energy function leads to bi-modal distribution

of bending angles

One clear consequence of non-convexity of the bending potential (Fig. 6), is that the corresponding distribution of polymer bend angles becomes bi-modal once the average bending angle $\bar{\theta}$ is within the convex hull region. A weak enough bending is always uniform, for as long as the average bend angle $\bar{\theta}$ is below $\theta_{a}$. As the average bend angle becomes just slightly larger than $\theta_{a}$, most of the segments are still bent weakly through $\theta_{a}$, but a small fraction becomes strongly bent through $\theta_{b}$. As the constraint forces the system to bend further, the fraction of the strongly bent segments increases linearly with $\bar{\theta}$, until, eventually all the segments are strongly bent through $\theta_{b}$. Beyond that point the system re-enters the uniform bending regime again. Here we confirm this expectation quantitatively, for a coarse-grained DNA model with one bead per bp. Specifically, we have analyzed angular distribution of different sized loops, from 6 to 600 bp long. The closed loops were created and simulated as discussed with the Main text. The loops were simulated at $T=300K$ to create 100,000 snapshots. The corresponding probability distribution of bend angles for each loop size is shown in Fig. 8.

Note that for a large, but finite $N\gg 1$, the predictions in Fig. 8 should be interpreted as applicable to the entire conformational ensemble of bent loops; in particular, the predicted linear dependence of the fractions of strongly and weakly bent fragments apply to the ensemble averages. For a given value of $\bar{\theta}$, between $\theta_{a}$ and $\theta_{b}$, one will observe a distribution of tightly bent fragments among the loops, a given loop can have none or more than one; bent conformations we observed in our MC simulations were qualitatively consistent with the above picture. In this work we did not pursue a detailed analysis of “structural” consequences of ECH model, in part because of an uncertainty associated with the value of $\theta_{b}$ appropriate for double-stranded DNA. We hope to revisit this issue in the future.

XI j-factor envelope functions and fitting to experiment

The difference between ECH and WLC models is encoded in the j-factor via

[TABLE]

where $\mathbfcal E_{loop}$ is the bending energy of forming a loop within each model, see Main Text. For quantitative comparison with experiment we follow the standard approach and augment the above formula with an oscillatory term that accounts for periodic variations in $J(L)$ due to torsional rigidity.

Specifically, the torsional dependence of the j-factor is represented by a cosine function, similar to that of the original work Shimada and Yamakawa (1984), to oscillate between the upper and lower envelope functions, with a period of $h=10$ bp, see Fig. 3 in the Main Text.

[TABLE]

where the “top”/“bottom” envelope curves have the functional form of the j-factor in Eq. 7,

[TABLE]

Here we are only interested in loops of length $L<160$ bp, since for $L>160$ bp ECH = WLC by construction for DNA. Thus, the curve representing ECH in Fig. 3 of the main text is Eq. 8, with the corresponding top/bottom envelope functions given by Eqs. 9.

For the case of WLC we replace $\mathbfcal E_{loop}(L)$ with the appropriate expression, so $J_{top/bot}$ used in Eq. 8 are now given by:

[TABLE]

The values of $k_{top},k_{bot}$ in the above equations are inferred from fitting Eq. 8 to experimental $J(L)$ points; since we have added an oscillatory part to $j(L)$ of Eq. 7, we need two parameters in the functional form of $J(L)$. For the fit we have chosen two experimental data points for loop lengths $L=101$ bp (the top envelope of $j(L)$) and $106$ bp (the bottom envelope). The reason we chose $J(L)$ that correspond to the largest $L$ available is because the resulting fit presents the most stringent test for the model in the tight bending regime where $L$ is small. Note that each curve in Fig. 3 of the Main Text, and Figs. 10(a), 10(b), has its own set of best fits values of $k_{top},k_{bot}$, which explains why the WLC and ECH curves do not coincide in the limit of large loops, where ECH approaches WLC.

The ECH envelope functions, Eqs. 9, have a minimum, which can be derived directly from the expression of j-factor (Eq. 7). The minimum occurs at

[TABLE]

Using the value of $L_{p}=150$ bp and $\theta_{a}=2.2^{\circ}$, the minimum of j-factor is found at $L\simeq 45$ bp. No such minimum exists in the WLC case in the range of $L$ of interest to us.

Note that ratios of J-factors at integer values of helical repeats can be predicted directly from Eq. 7, which does not contain the oscillatory components. The constant $k$ also cancels, which presents a convenient way to compare ECH theory directly to experiment, Table 1.

XII Robustness to model details

The key parameter of the ECH theory is the transition point, $\theta_{a}$ from the quadratic to the linear bending regime. As seen from Fig. 10, a good agreement with the experimental j-factors is achieved over a relatively broad range of $\theta_{a}$ values. Specifically, both $\theta_{a}=2.2^{\circ}$ inferred from experiment, and $\theta_{a}=1.5^{\circ}$ obtained from atomistic MD simulations (see Main Text) yield nearly the same agreement with experiment. In both cases, the counter-intuitive prediction of increased cyclization probability for very short loops hold.

Predictions of ECH framework are also robust to the precise value of the helical repeat parameter $h$ used in Eqs. 8 and 7 to account for the torsional stress created by loops with non-integer number of helical repeats. The use of $h=10.5$, more appropriate for DNA in solutionWang (1979), vs. $h=10$ corresponding to classical B-form, has little affect on the agreement of the predicted j-factors with the experiment.

XIII Origin of the non-convex bending at the atomic level

To illustrate, at the atomic level, the origin of the non-convex behavior of the backbone-backbone VDW energy, consider pairs of oxygen (O) atoms in the backbone. We choose these atoms because they have one of the lowest VDW energy minima out of all atom pairs in the DNA, and have most atom pairs within the effective short range of the interaction. The O-O VDW energy as a function of pairwise distance $r_{ij}$ has a potential well at $r_{min}\sim 3.3\AA{}$, Fig. 11. Further analysis of pairwise interatomic distances reveals that as the double helix bends, the geometry of the backbone deforms in a way that more oxygen atoms pairs fall into this well ($\sim 3.1-3.5$Å): beyond $\theta_{a}$ the accumulation of the attractive contributions begins to sharply lower the total interaction energy compared to the native, unbent state where all the oxygen atoms are significantly further apart ($>$4.0Å) than the well minimum. As the helix bends further, these O-O pairs eventually pass the LJ minimum, and begin to climb onto the repulsive wall, which increases the total VDW interaction energy, as seen in the corresponding figure of the main text.