A generalized theory of semiflexible polymers

DNA bending on length scales shorter than a persistence length plays an integral role in the translation of genetic information from DNA to cellular function. Quantitative experimental studies of these biological systems have led to a renewed interest in the polymer mechanics relevant for describing the conformational free energy of DNA bending induced by protein-DNA complexes. Recent experimental results from DNA cyclization studies have cast doubt on the applicability of the canonical semiflexible polymer theory, the wormlike chain (WLC) model, to DNA bending on biological length scales. This paper develops a theory of the chain statistics of a class of generalized semiflexible polymer models. Our focus is on the theoretical development of these models and the calculation of experimental observables. To illustrate our methods, we focus on a specific toy model of DNA bending. We show that the WLC model generically describes the long-length-scale chain statistics of semiflexible polymers, as predicted by the Renormalization Group. In particular, we show that either the WLC or our new model adequate describes force-extension, solution scattering, and long-contour-length cyclization experiments, regardless of the details of DNA bend elasticity. In contrast, experiments sensitive to short-length-scale chain behavior can in principle reveal dramatic departures from the linear elastic behavior assumed in the WLC model. We demonstrate this explicitly by showing that our toy model can reproduce the anomalously large short-contour-length cyclization J factors observed by Cloutier and Widom. Finally, we discuss the applicability of these models to DNA chain statistics in the context of future experiments.

bending. These protein-DNA interactions are of central importance to cellular function on a microscopic scale, from chromosomal DNA packaging, to transcription, and gene regulation, to viral packaging [1]. Protein-DNA interactions typically induce short-length-scale DNA bending which couples the chemical and physical properties of DNA [2,3,4].
A particularly important and successful application of polymer statistics has been in the description of double stranded DNA (dsDNA) by the wormlike chain model (WLC). In the WLC model, DNA is modeled as a fluctuating, linearly-elastic rod.
This simple model has been remarkably successful in describing many aspects of DNA mechanics and the statistics of semiflexible polymers generally. In particular, WLC describes the extension of a single dsDNA molecule under an external force with impressive precision [5].
Despite the notable theoretical and experimental success of the wormlike chain model, recent DNA cyclization studies by Cloutier and Widom [6] have cast doubt on the validity of the WLC model for describing the cyclization of short-contourlength sequences of DNA. In still more recent cyclization studies, Vologodskii and coworkers claim that the WLC model does accurately describe the cyclization of short DNA sequences [7]. Nevertheless, as we will explain later a number of experiments do seem to point to a role for elastic breakdown in DNA mechanics.
With the current experimental situation still in flux, it seems imperative to reevaluate the WLC model theoretically. We wish to answer the questions: (i) How could such a simple theory hope to describe a complex molecule like DNA? (ii) More precisely, which classes of experiments would we expect to be successfully described by WLC model, and which might require a different theory? Do these experiments correspond to the known successes or the recently reported failures of the theory?
In other words, we are asking how much room do the classic tests of WLC model leave for generalization of this model, and how completely do these experiments test the WLC model? Finally, we must ask (iii) Would a breakdown of the WLC model have any biological significance?
The focus of this paper will be the theoretical analysis of these questions and the development and discussion of more general semiflexible polymer models. Although these ideas are widely applicable to polymers statistics in general, the focus of this paper will be exclusively the mechanics of DNA. We shall attempt a synthesis of the existing experimental knowledge to determine which aspects of DNA bending are probed by existing experiments. In particular, we determine which experiments are most sensitive to the DNA mechanics relevant for understanding biological systems.
In the remainder of this introduction, we shall quickly outline our answers to the questions posed above.

A. Scale dependence in statistical physics
First, to put the possible breakdown of the WLC model into perspective, it is helpful to consider the bending of macroscopic rods. To engineers in the mechanics community, whose work has been the study of macroscopic bending, the failure of a linear elastic model at high curvature is more pedestrian than remarkable.
The linear elastic theory is understood to apply only to the small deflection limit.
What is perhaps more remarkable to some is that a linear-elastic model describes a macromolecular polymer at all, let alone to the accuracy illustrated by forceextension measurements! To put the success of the WLC model into perspective, it is helpful to consider DNA mechanics from the viewpoint of the statistical mechanics of condensed matter systems. Many physical properties of complicated condensed matter systems have been described by a small set of theories described in terms of renormalizable operators [8]. Regardless of the complicated structure of the theory at short length scales, the Renormalization Group (RG) guarantees that the long-length-scale chain statistics will be described by an effective energy functional containing only a few terms. In fact, for semiflexible polymers, only one such "renormalizable" term exists with the right symmetries. As a consequence, all semiflexible polymers share generic long-length-scale behavior: that described by the WLC model. Physically, this loss of information about the microscopic details of molecular mechanics arises from the averaging effect of thermal fluctuations.
The RG world-view leads us to expect that experiments like measuring the forceextension relation of long DNA would reveal only generic behavior, insensitive to microscopic details of DNA elasticity. But, on short enough length scales, the underlying structure of the theory becomes important. Violations of the linear elastic theory, analogous to those observed in macroscopic bending, are therefore expected in experiments that probe the short-length-scale bending of DNA. Indeed, early AFM imaging experiments did see the onset of deviations from WLC expectations on short scales [9]. Cyclization experiments, like the ones in Refs. [6,7], hold the prospect of greater sensitivity to the high-curvature regime.

B. Summary of this paper
This paper develops the qualitative framework outlined above by introducing a generalization of the wormlike chain model. This class of models, introduced in Sect. II, generalizes the WLC by describing a semiflexible polymer by an arbitrary local bending energy function. Sects. II C-II D introduce an explicit toy model of DNA bending, the "Sub-Elastic Chain" model (SEC), motivated by imaging data on DNA adsorbed to mica [10], and by recent nanoscale force measurements [11]. Sect. II E illustrates a computational procedure for computing the tangent distribution function for arbitrary contour length in generalized theories. Sect. II H introduces the persistence length in generalized theories and shows that these theories converge to the WLC model at long contour length.
The remainder of the paper focuses on the spatial distribution of the polymer.
The spatial distribution is of particular importance for biological applications where the contribution of chain statistics to biological function can often be formulated in terms of an effective end-concentration, the Jacobson-Stockmayer factor (J factor).
Physically, this effective concentration is the probability density of the polymer having the correct configuration for binding to the binding site of a protein. Sect. III introduces a method for computing the spatial and tangent-spatial distributions of generalized semiflexible polymer models in terms of a framework developed by Spakowitz and Wang [12,13,14] and others [15]. Sect. III A explicitly computes the spatial distributions for both the SEC and WLC models for various contour lengths. We discuss the Renormalization Group applied to spatial distributions and show the predicted convergence of the SEC and WLC models at long contour length. Sects. III B and III C show that the force extension and the structure factor computed for general theories are nearly identical to the WLC model results, implying that these experiments do not probe the high-curvature chain statistics important for many biological processes. Sect. III D computes the cyclization J factor for generic theories. We show that the SEC model gives rise to the enhanced cyclization efficiency for short-contour-length sequences observed by Cloutier and Widom [6] while leaving the long-contour-length J factor identical to that predicted by the WLC model. Finally, we discuss the results of this paper in the context of the recent cyclization measurements of Vologodskii and coworkers [7] and recent measurements of the deflection force for short sequences of DNA [11].
C. Relation to other recent work Following Crick and Klug's initial suggestion that DNA might kink [16], many classical works investigated the structural implications of this conformational change at the single-basepair level (reviewed in [17]). Indeed, many known DNAprotein complexes do display kinks in the DNA backbone. In contrast, our focus here is on physical measurements of DNA mechanics on a mesoscopic, severalbasepair scale relevant for biological processes like DNA looping in vivo. As described in Sect. II C, both Yan and Marko and we previously formulated and solved "kinking" models, in which DNA is assumed to undergo a sudden loss of bend stiffness beyond a critical stress [18,19,20]. Other related models were also formulated and solved in Refs. [21] and [22]. Sucato et al. have also performed Monte Carlo simulations of kinkable chains, to obtain information about their structural and thermodynamic properties [23]. Unlike these prior articles, the present work explores the proposal that the breakdown of linear elasticity, when coarsegrained to the mesoscopic scale, is effectively less abrupt than in the kinking models. We suggest that such a model can reconcile the growing evidence for elastic breakdown with the generic absence of sharply kinked states when tightly-bent DNA is observed microscopically.

A. Local energy functions
In this paper, we discuss a class of generalized elastic models for the statistical mechanics of semiflexible, inextensible polymers. The theories we discuss will be applicable to the description of polymers on length scales longer than the scale of the molecular structure. Accordingly, we idealize a semiflexible polymer as a chain (or "rod") consisting of N discrete segments ("links"), each of length ℓ, joined by semiflexible hinges ("vertices;" see Fig. 1). The link length ℓ should be taken to be shorter than the scale of the experiment we wish to describe, for example, shorter than the total length of the DNA in a cyclization experiment.
We then introduce a coarse-grained free energy cost for each chain configuration where E j is the energy associated with vertex j. To make a connection with the continuum mechanics picture, it is convenient to write this vertex bending energy as an energy density ε: that is a function of the N tangent vectors { t 1 , . . . , t N } and the vertex number j.
The coarse-grained configurational free energy E is a combination of entropic and energetic parts, which depend on the underlying molecular structure of the polymer.
We will ignore the effects of excluded volume, since we will be principally interested in bending on length scales where self-interaction effects play a negligible role in describing the chain conformation. We will also not allow for long-range interactions (that is, longer than ℓ); for instance, we assume that the solution conditions fix an electrostatic screening length smaller than ℓ.
To focus our attention on the novel effects of the hypothesis of elastic breakdown, we will restrict Eq. 1 to a subclass of models by making some assumptions about the form of the free energy density ε. Although this subclass is not a fully realistic depiction of known properties of DNA, it does have the virtue of being analytically The deflection angle between links i and i + 1 is θi.
tractable. Features of real DNA neglected in our models can be introduced in more numerical approaches once the phenomena we study are appreciated.
First we shall assume that the free energy density does not explicitly depend on the position j (the chain is homogeneous). Strictly, this is not the case for DNA since both the helical pitch and the sequence dependence spoil homogeneity [2,6].
We will study the mechanics of DNA on length scales longer than the helical repeat (3.6 nm), where helical effects approximately average to zero. Sequence dependence is a more serious omission [6], but we make this approximation in order to get analytical formulas. Having agreed to neglect helical effects, it is reasonable to add the assumption that the theory is rotationally invariant (the bending stiffness is isotropic). Last, we assume that the energy density involves only the first derivative of t, the discretized version of curvature: By rotational invariance, the energy density is a function of the magnitude of the curvature only: The most questionable of the assumptions above is the dropping of higher derivatives, which we will call "locality." For example, the role of nonlinear elastic-energy terms such as κ 4 will be a central concern of this paper, but this term has the same dimension (length) −4 as a term like (∆ κ /∆s) 2 , which we drop. Our justification is that higher-derivative terms correspond physically to cooperative conformational changes along the polymer, and although there are hints of long-range cooperativity in DNA [24], still the phenomena addressed in this paper do not seem to require such behavior. Again, our restricted class of theories is an analytically tractable starting point for the study of the effects of nonlinear-elastic behavior in an entropy-dominated chain. We will return to this point again briefly in the Discussion (Sect. IV).
An important one-parameter family of polymer theories is described by an energy density that is quadratic in the curvature: The models described by Eq. 5 have a restoring torque −dε/dθ that is linear in the link deflection θ: they are linear-elasticity theories. Their continuum limits are called wormlike chain models. A WLC model is completely characterized by one number, the elastic bending modulus ξk B T .

B. Statistical mechanics
We define the statistical mechanics theory associated to an energy function ε(s) in the canonical way. The probability of a coarse-grained chain conformation is given by the Boltzmann law: where Z is the partition function, determined by normalization, and we have defined the effective-coarse-grained free energy E in units of k B T .
The link length ℓ plays two key roles in Eq. 6. First, it defines the coarse-grained configuration space: a chain of physical length L tot consists of N = L tot /ℓ links, each with its orientation variable t j . Second, ℓ enters Eq. 4 explicitly. The physical meaning of ℓ is not obvious, however-it does not correspond to any crystallographic length in the DNA structure, for example the basepair rise of 0.34 nm. In fact, strictly speaking ℓ is not a parameter of the theory at all, because two different values of ℓ can give rise to theories with identical predictions, if the two theories' energy functions ε are suitably adjusted. Instead, ℓ is needed to give meaning to the energy function ε. The adjustment needed to maintain a fixed theory as ℓ is changed is called the "renormalization group flow" of ε [8].
It may seem tempting to eliminate ℓ from the theory by attempting a continuum limit, ℓ → 0, and indeed continuum mechanics does just this. When fluctuations are important, however, the continuum limit can discard some legitimate physics, and so must be treated with caution. In fact, we will argue that in the polymer context, the continuum limit leads only to a subset of models corresponding to the WLC (Eq. 5), because there is only one renormalizable term in the energy density with the right symmetries to meet our assumptions.
The fact that more general energy densities always have continuum limits describable by WLC models does not mean, however, that the WLC exhausts the physically legitimate and interesting models for stiff polymers. After all, we cannot expect continuum elasticity to remain valid on molecular length scales. Rather, this observation only implies that models more general than WLC must be defined with respect to some finite length scale ℓ.
The assumptions we made in Sect. II A imply that the partition function for an unconstrained N link chain decouples into N − 1 factors of the single-vertex partition function: where t i+1 is the out-going tangent of link vertex i (see Fig. 1). In the expression above, d t i+1 denotes an element of the d − 1-dimensional sphere of unit vectors in d = 2 or 3 dimensions; both two and three dimensions are of experimental interest.
We now introduce the fundamental tangent distribution function. The tangent distribution function is the conditional probability density for the final tangent, given an initial tangent. The fundamental tangent distribution function is the distribution function over just one link, length ℓ, and is related to the vertex energy by the Boltzmann Distribution (Eq. 6) where t i and t i+1 are the initial and final tangent respectively and the deflection angle at vertex i is given by cos θ i = t i · t i+1 . The chain statistics of the theory are completely determined by the fundamental tangent distribution.

C. Sub-Elastic Chain Model
We now introduce an explicit toy model for DNA bending that differs dramatically from the WLC model. Although the symbolic results below can be applied to the analysis of any of the semiflexible polymer models specified by Eq. 4, we will illustrate the method by using the toy model to compute experimental observables like force extension, the cyclization J factor, etc. for an explicit generalized theory.
The model we will study has a bending energy that is softer for high curvature than the WLC model. Nevertheless, it reproduces the successful long-length-scale predictions of WLC.
We have already described one such model in an earlier paper [19] and a similar model was also proposed by Yan and Marko [18,20]. In both cases, the high curvature softening was introduced by allowing kinking, or curvature localization: beyond a critical strain, the DNA's resistance to bending was supposed to fall suddenly to zero, or some small value. Although these kinking models reproduced the two desired features mentioned above, they predicted that highly curved DNA should be generically kinked [19]. However, atomic force microscopy (AFM) imaging of small minicircles generically shows them as round (although kinking can be induced in unusual ionic conditions) [25,26]. Moreover, tightly looped DNA shows enhanced sensitivity to DNAse digestion that is not concentrated on a single kink point, but rather is spread throughout the loop [27]. Finally, recent molecular-dynamics simulations of DNA minicircles show the spontaneous formation of sharp bends without strand separation [28]. For these reasons, this paper will explore a class of models with nonlinear DNA elasticity but without the catastrophic breakdown characteristic of kinking.
The bending energy functions we will study comes from the observation, well known in continuum mechanics, that a rod bending energy density that is nonconvex in curvature induces kinking when the rod is strongly bent [29]. To avoid kinking, we must therefore require that our effective bending energy density be everywhere a convex function of curvature, at least on length scales observable via electron microscopy (EM) or AFM imaging.
A simple choice of bending energy function that meets all the conditions mentioned, but is radically different from the WLC model, is: which defines a family of models parameterized by A and ℓ that we call "sub-elastic chain" (SEC) models. We will show that taking A = 5.3 and ℓ = 5 nm gives rise to a model with the persistence length ξ = 53 nm needed to describe the long-scale behavior of DNA in moderate-salt solution [30]. As a final motivation, AFM studies of the tangent-tangent correlation of DNA adsorbed to mica appear to fit a bending energy of roughly this functional form [10].
The SECmodel illustrates several of the points we wish to make. In particular, it is clear that high curvature, the bending stiffness in this SECmodel is softer than the corresponding WLC model at high deflection (we shall show that the persistence length of, and the energy is nowhere non-convex. Nevertheless, we will see that the SEC model reproduces essentially the same behavior as WLC for those aspects of DNA mechanics that have been well tested.
We emphasize that the SEC model defined by Eq. 9 is a toy model, and not intended as a realistic, accurate representation of DNA. In particular, the nonanalytic behavior of Eq. 9 at κ = 0 is not meant to be taken literally. Instead, it illustrates our calculational method, and our larger point that the classic DNAmechanics tests underdetermine even the coarse-grained effective DNA mechanics on biologically relevant length scales.

D. Measurements of the short-length-scale bending energy
The force required to tightly bend short sections of DNA has recently been directly measured by Liphardt and coworkers [11] via a fluorescence resonance energy transfer (FRET) force sensor. In this experiment, a sequence of DNA 9.18 nm in length is tightly bent by a linking sequence of single-stranded DNA as illustrated in Fig. 2. This contour length is represented in our theory by two links (ℓ = 5 nm) and a single vertex. The deflection angle is roughly 2π/3. It is straightforward to The deflection angle θ is roughly 2π/3 since the single-stranded DNA is roughly the same length as the link length.
estimate the deflection force in both the discrete WLC and SEC models: (In this estimate, we have used ℓ = 4.6 nm, half the contour length of the dsDNA.) The experimentally measured force, 6 ± 5 pN, is approximately equal to that predicted by the SECmodel but is more than a factor of two smaller than that predicted by the elastic rod model (WLC). These experiments again indicate that the WLC model fails to describe the high-curvature bending of short sequences of DNA. At least at this deflection angle, the SEC model approximately predicts the deflection force. Note that if the kinking model of Refs. [18,19] literally described short sequences of DNA, this force would be zero, contrary to the experiment-another motivation for our introduction of generalized elasticity models.

E. The propagator and composition
The locality assumption in the definition of the bending energy implies that each vertex bends independently. The fundamental tangent distribution function is the conditional probability of a final tangent, given an initial tangent for a single vertex. Computing the tangent distribution functions for chains of several links is therefore straightforward. These conditional probabilities are simply the product of conditional probabilities for single vertices, summed over the orientations of the intermediate tangents [20] where we have written the N link tangent distribution function as a function of the arc length, N ℓ. This notation is needlessly cumbersome. It is therefore convenient to introduce the propagation operator (or transfer matrix [20]) where | and | is the canonical bra ket notation of statistical mechanics (or quantum mechanics) [31]. These states are a continuum basis: where δ is the Dirac delta function on the space of unit tangent vectors.
The propagation operator, G, applied on a state gives the state (probability distribution) after one additional link. This property is called composition and is a direct result of the locality discussed above. We can now rewrite Eq. 12 more concisely where the weighted sum, or path integral, over all intermediate configurations is now implicit. By changing the basis in the next section, we shall show that this expression is also a convenient computational tool for understanding general theories [20].

F. Symmetry considerations
The tangent basis we have exploited to write the tangent distribution function is not particularly convenient computationally since the operator is not expressed in its eignenbasis in which it is diagonal. To find an eigenbasis for this operator, we exploit the rigid body rotational invariance of the tangent distribution function. In D dimensions, the rigid-body-rotational invariance of the model implies that the propagator commutes with the generators of rotation where L ij = −L ji are the generators of rotation in the ij plane. The propagator therefore also commutes with the Casimir operator, which in Quantum Mechanics would correspond to the total angular momentum: Since L 2 and G commute, they share an eigenbasis [31]. The angular momentum states span the tangent space and are eigenvalues of L 2 : where l is the total angular quantum number and we write the other angular quantum numbers collectively as m. The propagator can therefore be expanded in this where the g l are coefficients that depend only on the quantum number l but not on m. Eq. 19 is the desired diagonalization of the propagator G.
Explicitly, in two dimensions, it is convenient to use the eigenfunctions [31] for integer m and l ≡ |m|. Note that the quantum number m is sufficient to describe the state but we have introduced a second quantum number, l, which is invariant under a generalized notion of rotational invariance in two dimensions, including the discrete transformation θ → −θ (parity inversion).
In three dimensions, it is convenient to use the eigenfunctions [31] t where the Y lm are the spherical harmonics. In this case m is the eigenvalue of the z component of the angular momentum operator L 12 .
The orthonormality of the basis implies that the g l are uniquely determined and can be found in the usual way (Appendix A Eq. A1 and Eq. A2). It is now straightforward to perform the N+1 link path integral of Eq. 12 since the propagation operator is diagonal.
We return now to the SEC model proposed in Sect. II C. The N link tangent distribution function for the SECmodel is shown in Fig. 3. This figure explicitly illustrates the scale dependence of statistical mechanics theories. For short-contourlength chains, the WLC and SECtheories make dramatically different predictions, but as the contour length of the chain increases, the differences between the distribution functions of the two theories decrease until at long contour length, the theories are essentially indistinguishable. This is the essence of the renormalization group: at short length scales, the mechanics of the chain can be extremely com-

G. Contour length continuation
Since we will frequently be interested in the properties of the polymer on length scales much longer than the fundamental link length ℓ, it is useful to introduce a

Hamiltonian operator defined by
H is also diagonal in the angular momentum representation with eigenvalues h l = −ℓ −1 log g l . We call this operator the Hamiltonian operator because in the WLC model, the statistical mechanics of the chain corresponds to a quantum particle on a D − 1 sphere. The tangent distribution function is equal to the quantum propagator where time has been continued to imaginary arc length. The operator H is equal to the Hamiltonian of the corresponding quantum mechanical particle system.
We can rewrite the N link as propagator where N ℓ is the contour length corresponding to N links. The advantage of this reformulation of the distribution function is that it introduces a natural extension to fractional numbers of links by replacing N ℓ by the contour length L defined for all positive real numbers: although rigorously, it is understood that this function is only defined for contour lengths equal to an integral number of links.
H. The meaning of persistence length, and the stiff-polymer limit What is the meaning of persistence length in general models like the ones we have described? Persistence length describes the length scale on which the polymer maintains its tangent orientation. For the WLC model in D dimensions, the tangent persistence is where ξ is the persistence length, which also appears in the energy density as the bending modulus in Eq. 5. In general models, the tangent persistence (Eq. 26) has the same functional form but ξ no longer corresponds to a bending modulus. We shall therefore simply use Eq. 26 to define the persistence length of general models.
The tangent persistence corresponds to the l = 1 mode of the propagator. (In the quantum mechanical correspondence, t is a vector and creates a state of spin one.) Comparing Eqs. 25 and 26, the persistence length is where h 1 is the l = 1 eigenvalue of the Hamiltonian. Note that we have explicitly written a subscript D to denote the dependence on dimension. In the WLC model, ξ is independent of dimension, but in more general models this is not the case.
The persistence length also controls the long-length characteristics of the polymer. The mean-squared end-to-end distance can be written in terms of the tangent persistence Since Eq. 26 applies to both the WLC model and general models, the dependence of the mean squared end-to-end distance on persistence length and contour length is identical to the WLC model. The same is true for radius of gyration, which can also be written in terms of an integration of Eq. 26. It is also well known that the long-contour-length spatial distribution of semiflexible polymers is described by the Gaussian Chain model [32]. The width of the Gaussian distribution is determined by the mean squared end-to-end distance; the relation between the Kuhn length and the persistence length is therefore the same for our general models as for the WLC model.
We can immediately exploit Eq. 27 and Eq. A2 to analyze the SECmodel. The persistence length, computed for the SECmodel in three dimensions is 53 nm which matches solution measurements.
Stiff-polymer limit: We now examine the tangent distribution function in the stiff polymer limit and show that the WLC model is universal at long contour length as predicted by the Renormalization Group. Our explicit computations of the SECtangent distribution function in Sect. II E have already provided one explicit example of this behavior, but we address this question generally in this section.
By definition the stiff polymer limit implies that the fundamental tangent distribution function, g, is narrowly distributed around zero deflection. We will exploit this fact by expanding the basis functions in the deflection angle and computing the eigenvalues of the propagator (Eq. 19) to lowest order in the deflection angle.
In dimension D, this calculation, although straightforward, requires some technical mathematics. We therefore relegate the details of this calculation to the appendix, Sect. B, and present only the results.
The propagator in the stiff polymer limit is Eq. B13 where ξ is the persistence length defined by the l = 1 eigenvalue of the Hamiltonian operator (Eq. 27). Note that the L 2 term is understood to be small for small values of the angular quantum number l since, in the stiff polymer limit, the link length ℓ is much shorter that the persistence length ξ. The corrections are order L 4 θ 4 and scale as l 4 for large l. Clearly this approximation holds only for small angular quantum number l. It is convenient to compute the Hamiltonian operator which is identical to the WLC Hamiltonian to lowest order in the deflection angle.
Again, the correction scales like l 4 which implies that this relation holds only for small angular quantum numbers.
The correspondence between the Hamiltonian operators for general models and the WLC model for small angular quantum numbers implies that the long-contourlength behavior of the polymer is universal and determined by the persistence length alone. This correspondence is shown explicitly for the SECand WLC theories in The details of the short-length-scale bending energy affect only the large l eigenvalues of the Hamiltonian operator and are therefore irrelevant at long contour length, as predicted by the renormalization group.
Although we have yet to compute the spatial distribution function, we have explicitly shown that measurements that are only sensitive to the long-length-scale

III. THE SPATIAL DISTRIBUTION
For most applications, it is the spatial distribution of the polymer rather than the tangent distribution function which is of phenomenological interest. From solution scattering to force-extension, cyclization to looping, the spatial distribution function is directly observable. In this section, we shall develop a near exact formalism for computing the spatial distribution function. Our focus will be exclusively three dimensions but computations in other dimensions are a simple extension of the methods discussed here.
The tangent-spatial distribution function is the probability density of end displacement X and final tangent t f given an initial tangent t i for an arc length L chain. It is convenient to write the tangent spatial distribution in terms of the spatial delta function [33] G( X; t f , t i ; where we have written the distribution function in the continuum limit. We shall reintroduce an operational definition of this continuum limit in a moment.
To compute the tangent-spatial distribution function, we introduce an operatorvalued spatial distribution function [13]: which allows us to keep the tangents implicit in our expressions. We shall call this operator the spatial propagator since it obeys the composition property of Green Functions: where ⊗ denotes the spatial convolution.
As before, it will be convenient to work in the angular momentum basis with the matrix elements since this basis diagonalizes the Hamiltonian (although not the spatial propagator).
Finding the spatial propagator reduces to the ability to explicitly compute all the We shall be able to derive exact expressions for the Fourier-Laplace Transform of the spatial propagator in the continuum theory in terms of the transformed matrix elements (Eq. 36). We adopt the Fourier Transform convention and the Laplace transform conventioñ The derivation of the transformed matrix elements exploits the same techniques used recently by Spakowitz and Wang [12,13,14] to derive exact results for the WLC model. The extension of these results to the generalized theories considered here is straightforward. We shall therefore include only a brief derivation in Appendix C although we discuss the results in the main text.
It is important to note at this stage that the results derived for the spatial distribtution, although derived in a method analogous to that exploited in Ref. [13], will not be exact solutions to generalized discrete-link models. Rather, the results quoted here are exact-solutions to the analytically continued theories defined by Eq. 25. That is, we have assumed a formulation of the discrete-link theories that defines the tangent distribution function for all L > 0, although formally this distribution function is defined only for contour lengths equal to an integral number of links. For semiflexible polymers longer than a few links, this is an excellent approximation. (For instance, the discrete and continued theories are later compared in Fig. 5.) We have therefore called this solution "near-exact" in the text.

A. The spatial distribution function
In force-extension and solution scattering experiments the tangents of the polymer are not directly probed by experiment; it is only the spatial distribution function rather than the tangent-spatial distribution function which is observed. We shall therefore introduce the spatial distribution function, K( X; L), which is defined as the probability density that a contour length L polymer has end displacement X. The spatial distribution function is the tangent-spatial distribution function summed over the final tangent and averaged over the initial tangent: where the last equality expresses the spatial distribution function in terms of a matrix element of the spatial propagator.
The Fourier-Laplace transform of this matrix element, a continued fraction, is computed in Eq. C15. The explicit expression for the transformed spatial distribution function isK ( k; p) = 1 where the h l are the eigenvalues of the Hamiltonian operator, Eq. 23, and the B n coefficients are defined as This expression is identical to that derived for the WLC model [12], except that the eigenvalues of the Hamiltonian operator, h l , are those for a generic theory rather than the WLC eigenvalues. Otherwise the expression is unchanged.
The spatial distribution functions for the WLC and SEC models are plotted in for small wave number, implying that the higher-angular-momentum states have successively less relevance at long length scales.
It is also instructive to consider the long-length-scale limit of the spatial distribution function since we know that this limit is described by the gaussian chain model. The long-length-scale limit corresponds to the limit of small k and contour dual number p. In this limit, the transformed spatial distribution function is which is just a Gaussian distribution with a Kuhn length of twice the persistence length (Eq. 26) as we have already argued from computations of the mean-squared end-to-end distance and has also been shown schematically for the SEC model in Fig. 5.

B. Force-extension
The force-extension of single polymer molecules has long been the subject of experimental interest [5,34]. force is induced by the reduction in the number of micro configurations available to the chain as the extension is increased.
The successful comparison of WLC to single-molecule force-extension data has been described as the strictest test of the WLC model [5]. But how do other semiflexible polymer models compare? Can these models also reproduce the precise fit to experiment? To answer these questions, we next compute the force-extension for general models and explicitly compare the extension in the SEC and WLC models (Fig. 6).
The partition function for a polymer under a constant external tension is related to the Fourier transform of the spatial distribution function via an analytic continuation of the wave number: where f is the external force or tension. The extension or mean end distance is computed in the usual way: The force-extension for the SECand WLC models are compared in Fig. 6. The numerical technique applied in this computation is described in Appendix D.
Despite the drastically different bending energy of the SEC model on short length scales, thermal fluctuations disguise these differences and give rise to an extension almost identical to the WLC model. In retrospect, these results are hardly surprising. The theories are identical at small extension due to the renormalization group and at large extension due to inextensibility. Although, in principle, the high force limit is mathematically equivalent to probing short length scales-they are related by analytic continuation-these differences are not large enough to be experimentally observable. Physically, the rare high curvature bending regime, where the difference between the models is most pronounced, is further suppressed by the application of tension. For the study of DNA mechanics, force-extension measurements do probe the persistence length and the inextensibility of DNA, but these experiments do not effectively probe DNA elasticity on the length scales of interest for many biological processes.

C. Structure factor
Another experimental observable used to characterize polymers is the structure factor, measured by static light scattering, small-angle X-Ray scattering, and neutron scattering experiments. Measurements of the structure factor can probe the polymer configuration on a wide range of length scales. Symbolically the structure factor is where X(s) is the position of the polymer at arc length s and we have included an extra factor of the contour length in the denominator to make the structure factor dimensionless [12]. At high wave number, the structure factor is sensitive to short-length-scale physics, whereas the contour length and radius of gyration are probed by low wave number. The structure factor can be rewritten in terms of the transformed spatial distribution function where L −1 is the inverse Laplace transform which can be computed numerically.
As mentioned above, the leading-order contributions at small wave number are determined by the polymer length and the radius of gyration where we have temporarily restored the length dimension of S. At large k, both WLC and SECare straight, which gives an asymptotic limit for large wave number The structure factor is compared for the SEC and WLC models in Fig. 7. The numerical technique applied in this computation is described in Appendix D. Again we find that the two theories make nearly identical predictions. The reasoning is again similar to that explained for force-extension. The two theories make dramatically different predictions for rare, highly-bent configurations but the structure factor is dominated not by these rare high curvature configurations but by typical thermal bending. We therefore find that the structure factor, like force extension, does not effectively probe the high-curvature statistics of the polymer.

D. Cyclization
The biochemical process of DNA cyclization is not in itself a process of particular biological importance 1 but cyclization experiments do provide a controlled, bulk experimental method for probing the probability of rare, highly-curved DNA configurations [6,7,37]. polymers [4]. The ratio of the cyclization equilibrium constant (K C ) to the dimerization equilibrium constant (K D ) is called the Jacobson-Stockmayer factor [38] or J factor and is proportional to the tangent-spatial distribution function of the polymer [4,40] J ≡ K C /K D = 4πG(0; t, t; L) = tr G(0; L), where G is the tangent-spatial distribution function for end-to-end displacement 0 and aligned end tangents, for a contour length L polymer. The J factor can also be written as the trace of the spatial propagator. (The matrix elements of the spatial propagator are written explicitly in Appendix C.) Physically, the J factor is proportional to the concentration of one end at the other with the correct (aligned) orientation for hybridization.
Our analysis neglects the condition that DNA twist must also be aligned, which requires the use of models including the twist degree of freedom. This additional constraint modulates the J factor with a 10.5 bp period equal to the helical repeat.
Our interest here is in the value of the J factor averaged over a helical repeat for which the effects of twist can be roughly ignored [33]. The renormalization group predicts that the SECmodel will be identical to WLC for longcontour-length sequences. But, for sequences shorter than two persistence lengths ( 200 bp), the short-contour-length chain statistics become important and the SECJ factor diverges from the WLC prediction. In fact, for 94 bp sequences, the SEC J factor is three orders of magnitude larger than that predicted by the WLC model, roughly matching the J factors measured by Cloutier and Widom [6,37]  Our results predict that a shortcontour-length anomaly in the J factor is generic for sufficiently short sequences, but the contour lengths at which the WLC model breaks down is model dependent.
applications and (ii) cyclization experiments are extremely sensitive to the differences between models at short contour length.

E. Beyond the J factor
The J factor is not the only effective concentration of interest. DNA looping is integral to the function of many gene regulatory proteins. The affinity of these proteins for DNA, and therefore their function, depends sensitively on the looping free energy, or equivalently the effective concentration of the looped DNA. (For instance, see Refs. [43], [44], and [45].) Once the geometry of the loop is known-the displacement of the binding sites ( X) and the orientation of the bound DNA ( t and t ′ )-both the SECand WLC models make predictions for the effective concentration: These statistical mechanics predictions can then be directly compared with quantitative measurements of gene expression [44,45] and in vitro experiments [46].
More general cyclization measurements may also be performed. For instance, cyclizing sequences with two short single-stranded gaps could be used to probe short-length-scale DNA mechanics. The short single-stranded sequences are very flexible and can be approximated as free hinges. This technique could be exploited to directly measure the spatial distribution function shown in Fig. 5 for very short sequences of DNA.

A. The SEC
Sects. II C-II D introduced the SEC as a toy model for DNA bending, motivated by several physical measurements on DNA. We proceeded to show that this simple model exhibited the long-length-scale chain statistics of the WLC model, despite dramatically increasing the predicted probability of high-curvature configurations.
In particular, we showed that the SEC model yields a cyclization J factor in agreement with the measurements of Cloutier and Widom [6]. More generally, we argued that this type deviation from WLC behavior is generic in semiflexible chain models.
These putative deviations of DNA chain statistics from the wormlike chain model at short contour length are quite relevant for structural biology, where the typical radius of curvature induced by DNA binding proteins is on order nanometers or tens of nanometers, not persistence lengths. For example, the radius of curvature of DNA bound in a nucleosome complex is roughly 6 nm. The structure of this complex shows sharp bends, but no sign of melting, consistent with our SEC model [47]. Similarly, DNA looped by a gene regulatory protein is typically bent on short length scales [43]. If DNA is described by the SECmodel, these tightly-bent DNAprotein complexes are orders of magnitude more stable than predicted by the WLC model. A quantitative understanding of biological DNA bending therefore awaits a consistent model of short-length-scale DNA bending.
Unfortunately, precise quantitative tests of short-length-scale DNA bending are still in the future. Vologodskii and coworkers recently made measurements questioning the results of Cloutier and Widom [7]. Their measurements suggest that the J factor agrees with that predicted by the WLC model, at least down to a contour length of 100 bp. Widom and coworker then repeated their own measurements, however, and have confirmed their previous results [48]. Also, Sect. II D mentioned that Shroff et al. [11] [49]. The twist density needed to join a mis-phased DNA loop of under 100 bp exceeds this threshold, and indeed Cloutier and Widom have also shown that the twist-induced modulation of the cyclization J factor is smaller for short sequences than predicted by the HLWC model [37].
Thus, although the bending of DNA for small twist densities may be adequately described by the HWLC model, a generalized model of DNA, including elastic breakdown of both bend and twist stiffnesses, may be necessary to describe the chain statistics of short sequences of looped DNA that are not naturally in twist registry when bound. Such generalized models are in principle a straightforward extension of the theory presented in this paper and new exact results for the HWLC model recently derived by Spakowitz [14].

V. CONCLUSION
In this paper, we have explored a class of generalized semiflexible polymer models in which the bending energy density is an arbitrary function of curvature. To analyze the chain statistics of these models, we developed a formalism that is analogous to the techniques used for describing the WLC model. We demonstrated that the statistics of these general models coincide with those of the linear-elastic (WLC) model at long contour length, as predicted by the renormalization group. At short length scales, we show that the predictions of these models can be dramatically different from the WLC model. We computed near-exact expressions for the transformed spatial and tangent-spatial distribution functions with a method analogous to that recently exploited to find exact results for the WLC model. These generalized models provide an explicit example of a non-renormalizable model which is nearly exactly solvable. We exploited these general theoretical results to compute several important experimental observables: force-extension, the structure factor, and the cyclization J factor. We explicitly performed these computations for a toy model of DNA bending, the Sub-Elastic Chain (SEC) model. The predictions of this model are essentially indistinguishable from the WLC model for force-extension, solution scattering, and long-contour-length cyclization measurements, despite dramatic differences between the bending energies of the two models on short length scales. For short-contour-length cyclization experiments, general models generically predict large deviations from WLC behavior. In particular we computed the J factor for the SECmodel and showed that this model could account for the anomalously large cyclization J factor measured by Widom and Cloutier [6]. We expect these generalized models to be widely applicable for describing the high-curvature statistics of other semiflexible polymers.
where θ is defined as the angle away from the z axis: t(0) = e z . In three dimensions, the g l are where the P l are the Legendre Polynomials and cos θ = t · e z . We begin the derivation with the definition of the lth moment of the tangent distribution function expressed in terms of the propagator Eq. 19 where rigid-body-rotational invariance implies that g l is independent of m. We insert two complete sets of states into the tangent representation We can now replace the matrix element of the propagator with the fundamental tangent distribution function g( t; t ′ ) (Eq. 13). Remember that this function depends only on the relative deflection angle of the tangents. We therefore replace the integral over the second tangent with an integral over rotation matrices, R, and make the substitution t ′ ≡ R t: where we represent the change in measure symbolically and we have introduced the rotation operator [31] Our interest is in the case where the tangent distribution function is narrowly distributed. We shall therefore expand the rotation operator, D, with respect to the rotation angles which we shall assume are small. The rotation operator can be expanded in terms of these angles and the rotation generators [31] where the θ ij = −θ ji are the components of the rotation angle which multiply the generators of rotations in the ij plane.
To evaluate the integral over the rotation matrices, we must now choose a set of θ's which give a single cover of the tangent space. Since g( t; R t ) is independent of t, it is convenient to choose a coordinate system in which t is in the direction of the D axis. (We shall return to the unrotated frame before performing the integral over t.) In this new coordinate system, it is convenient to use the cover generated by the coordinates {θ Di } 1..D−1 while setting all other θ's to zero.
We denote the average taken with respect to the distribution function by .
Due to rigid-body-rotational invariance around the D axis, where is the total deflection angle.
The nonzero matrix elements can be put in a coordinate invariant form lm| e D e D | L Di L Di |lm = lm| e D e D | L 2 |lm (B10) since the added terms in the Casimir operator, L 2 , are zero on | e D . We can now go back to the unrotated coordinate system by setting e D = t.
After integrating over the complete set of tangent vectors, the resulting moment is Since this expression is only correct to O(θ 4 ), it is convenient to replace 1 2 θ 2 with 1 − cos θ. We can now use the definition of the persistence length given in Eq. 26 to eliminate the dependence on cos θ : Finally, we reconstruct the propagator from its moments which completes the derivation. This result is discussed in Sect. II H.

APPENDIX C: THE TRANSFORMED SPATIAL PROPAGATOR
To derive closed form expressions for the spatial propagator, we Fourier Transform the spatial propagator over the relative displacement, X. In particular, we consider the Fourier Transform of Eq. 35 since in Fourier space, the spatial convolutions are simply products: We choose a coordinate system where k is in the z direction.
We now wish to use this composition property of the spatial propagator to write a differential equation for G. We therefore consider G for a differential arc length dL and then expand the Fourier Transform of Eq. 33 for arc length dL: where I is the identity operator and A ≡ H + ik cos θ where θ takes its canonical meaning in spherical coordinates: cos Θ = t ·ẑ. Substituting this expression into Eq. C1, we can write a differential equation forG: It is now convenient to make a Laplace transform from arc length L to its conjugate variable p. After solving for the propagation operator, we have an operator equation for the Laplace-Fourier Transform of the spatial propagator: but this expression is not explicit since it is written in terms of the inverse of an infinite dimensional operator.
We can express cos Θ in the angular momentum basis. It is most convenient to define a set of ladder operators: where the ladder operators are defined by and the A l,l+1,m are: A l,l+1,m = A l+1,l,m = (l − m + 1) (l + m + 1) (2l + 1) (2L + 3) .
The ladder operators have the property that they increase (decrease) the total momentum quantum number of a state by plus (minus) one.
Next, we obtain explicit expressions for the matrix elements of the transformed spatial propagator. The Hamiltonian is diagonal in the angular representation, so it is convenient to factor the spatial propagator (Eq. C4) into diagonal and nondiagonal factors: and expand it in a power series As a first step, we will compute a diagonal matrix element: Computing these matrix elements is achieved by grouping the infinite set of terms in Eq. C10 into sub sets which can be summed exactly [12].
We introduceG + l ′ ml ′ m which is the matrix element of a subset of the terms in Eq. C10, in which there are only transitions to states with total momentum l = l ′ or greater [12]. This matrix element can be defined recursively since only transitions to adjacent states are possible. The matrix element is the sum over n of the matrix elements with n transitions to and from the l ≥ l ′ + 1 states, which can be written in terms ofG + l ′ +1,m,l ′ +1,m . The terms of this matrix element, a geometric series, can be summed exactly [13]: This sum is pictured schematically in Fig. 9.
Similarly, we defineG − l ′ ml ′ m which is the matrix element of the propagation operator which allows transitions to states with total momentum l = l ′ or less: In terms of G ± we can now define the matrix element without transition restrictions by grouping the transitions into sets that do not cross l = l ′ . These sets can be written in terms of the matrix elements of G ± and then summed in a geometric series [13]: The diagonal matrix element computed above is sufficient for describing many observables of phenomenological interest. Note that the only difference between this expression and the WLC expression [13] is that the eigenvalues of the Hamiltonian operator have changed. This matrix element is represented by the line with ellipses, representing all transitions to states with higher l. (b) G + lmlm can be defined recursively in terms ofG + l+1,m,l+1,m . The definition ofG − lmlm is analogous, but it is the sum of all diagrams with transitions to states with total angular quantum number l and smaller.
For some applications we will want completely general matrix elementsG lml ′ m ′ .
We can again define these general matrix elements in terms of the recursive definitions of G ± . Again, the trick to summing the terms is grouping them. In this general case, there are many equivalent ways of achieving this grouping. See Fig. 10 for an explanation of the set grouping. The matrix element can be written [13] G l+n,m,l,m ′ =G l,m,l+n,m ′ = δ m−m ′G lmlm n q=1 −ikA l+q−1,l+q,mG + l+q,m,l+q,m .

(C16)
We have now explicitly solved for spatial propagator having written expressions for all the matrix elements.

APPENDIX D: THE COMPUTATION OF SPATIAL DISTRIBUTIONS
The previous section discussed near-exact expressions for the Fourier-Laplace transformed spatial and tangent-spatial distribution functions. Exact closed-form expressions for these functions are unknown and we must invert the transforms numerically to compute the distribution functions.

Force-extension and the structure factor
The computations of force extension and the structure factor require only a single numerical inverse Laplace Transform. We cut off the continued fraction at l = 10 FIG. 10: General matrix elements. A diagram of the sum for the matrix element G lml+nm =G l+nmlm . To compute the matrix element, we group the terms by the location of the first steps from l + n to l + n − 1 and from l + n − 1 to l + n − 2 etc. In the diagram, these steps are represented by the vertical lines. We use the G + operator to sum over all possible diagrams with upward transitions between these steps. These upward transitions are represented by the ellipses. We multiply by the transition matrix element for each of the vertical lines. After we reach l for the first time we allow all transitions up or down. This enumeration counts each contributing diagram once but this recipe is not unique. and then used the InverseLaplaceTransform function in Mathematica.

The spatial distribution and the J factor
For computations of the spatial distribution function and the J factor, we exploited two different numerical techniques: numerical transform inversion and Monte Carlo. For contour lengths of a persistence length and above, it is convenient to directly invert the transforms numerically by truncated the continued fraction in the transformed propagator (Eq. C15). Typically we used l = 15 as the cutoff although in some cases higher l values were used for short contour lengths.
In the inverse transform technique, both numerical Laplace and Fourier Transform inversions must be computed. We have used two different implimentations for these computations. (i) In Mathematica, we used the InverseLaplaceTransform function. We then integrated numerically (using an explicit sum) to invert the Fourier transform. We found that the built-in numerical integration in Mathematica was too slow for practical use. For contour lengths on order a persistence length and shorter, inverting the transformed expressions is impractical. The continued fraction in increasing momentum is essentially an expansion around weak end-tangent correlation. For contour lengths shorter than a persistence length, a larger l cutoff is required, significantly slowing the numerical inversions. In addition, the numerical integration over the wave number becomes impractical since the numerical integrations must be ex-tended to very a large cutoff momentum. These convergence issues are not unique to the continued fraction approach. For example, the transfer matrix approach is plagued by similar shortcomings, requiring difficult numerical work at short contour length [18].
We therefore used a much simpler, although less elegant, solution in the form of direct Monte Carlo integrations. Monte Carlo integration in the short-contourlength regime (i) is numerically more efficient than direct inversion, (ii) requires very minimal implementation, and (iii) serves as a useful check of our theoretical results. These checks appear few places explicitly in the paper since the agreement between these two methods is excellent and the focus of this paper is physics rather than numerical computations. The theoretical curve for the cyclization J factor