Unraveling looping efficiency of stochastic Cosserat polymers

Understanding looping probabilities, including the particular case of ring closure or cyclization, of fluctuating polymers (e.g., DNA) is important in many applications in molecular biology and chemistry. In a continuum limit the configuration of a polymer is a curve in the group SE(3) of rigid body displacements, whose energy can be modeled via the Cosserat theory of elastic rods. Cosserat rods are a more detailed version of the classic wormlike-chain (WLC) model, which we show to be more appropriate in short-length scale, or stiff, regimes, where the contributions of extension and shear deformations are not negligible and lead to noteworthy high values for the cyclization probabilities (or J-factors). We therefore observe that the Cosserat framework is a candidate for gaining a better understanding of the enhanced cyclization of short DNA molecules reported in various experiments, which is not satisfactorily explained by WLC-type models. Characterizing the stochastic fluctuations about minimizers of the energy by means of Laplace expansions in a (real) path integral formulation, we develop efficient analytical approximations for the two cases of full looping, in which both end-to-end relative translation and rotation are prescribed, and of marginal looping probabilities, where only end-to-end translation is prescribed. For isotropic Cosserat rods, certain looping boundary value problems admit nonisolated families of critical points of the energy due to an associated continuous symmetry. For the first time, taking inspiration from (imaginary) path integral techniques, a quantum mechanical probabilistic treatment of Goldstone modes in statistical rod mechanics sheds light on J-factor computations for isotropic rods in the semiclassical context. All the results are achieved exploiting appropriate Jacobi fields arising from Gaussian path integrals and show good agreement when compared with intense Monte Carlo simulations for the target examples.


I. INTRODUCTION
It is widely known that polymers involved in biological and chemical processes are anything but static objects. In fact, they are subject to stochastic forcing from the external environment that lead to complex conformational fluctuations. One of the fundamental phenomena which is understood to perform a variety of roles is polymer looping, occurring when two sites separated by several monomers, and therefore considered far from each other, come into proximity. A basic observation is that the interacting sites alone do not characterize the phenomenon of looping, but rather it is the whole polymeric chain that rearranges itself for this to occur. As a consequence, the length and mechanical properties of the chain, together with the thermodynamic surrounding conditions are finely tuning the likelihood of such events. There are many reasons to study this topic, which have led to a considerable literature. For instance, looping is involved in the regulation of gene expression by mediating the binding and unbinding of DNA to proteins [1][2][3], such as the classic * corazzagiulio24@gmail.com example of the Lac operon [4,5]. In addition, DNA packaging (chromatin formation) [6], replication and recombination [1,7] depend on the ability of the polymer to deform into loop configurations, as do other cellular processes. Proteins exhibit intrachain loops for organizing the folding of their polypeptide chains [8]; e.g., antibodies use loops to bind a wide variety of potential antigens [9]. When dealing with a closed loop, it is usually appropriate to refer to cyclization or ring closure. In this regard, the production of DNA minicircles is being investigated for their possible therapeutic applications [10]. Even in the context of nanotechnologies, ring closure studies have been performed for carbon nanotubes subject to thermal fluctuations [11] and wormlike micelles [12].
From the modeling point of view, it is appropriate to look back at some of the historical milestones that underpin our work. In 1949, Kratky and Porod [13] introduced the wormlike-chain (WLC) model for describing the conformations of stiff polymer chains. Soon after, the complete determination of the polymeric structure of DNA guided scientists towards the application of WLC-type models in the context of DNA statistical mechanics, allowing probabilistic predictions of relevant quantities of interest. Historically, the computations have been performed in terms of Fokker-Plank equations [14,15], but also exploiting the point of view of path integrals [16][17][18], a technique inherited from Wiener's work [19,20] and quantum mechanics [21]. These ideas were largely investigated by Yamakawa [22][23][24][25][26][27][28][29], who in particular considered the problem of computing ring-closure probabilities, now ubiquitous in molecular biology [30][31][32]. Nowadays, for a homogeneous chain, the exact statistical mechanical theory of both the WLC and the helical WLC (with twist) is known [33][34][35], and the topic has been rigorously phrased over the special Euclidean group SE(3) [36].
In parallel, in the early years of the 20th century, the Cosserat brothers Eugène and François formulated Kirchhoff's rod theory using what are now known as directors [37]. However, the difficulties arising from the generality of the model, which includes the WLC as a particular constrained case, hindered its application to stochastic chains. Only quite recently, targeting a more realistic description of DNA, the mentioned framework has been partially or fully exploited both within new analytical studies [38][39][40][41][42][43] and intense Monte Carlo (MC) simulations [44][45][46][47], the latter being only a partial solution because of time and cost.
In this article we aim to fill the gap between user-friendly but simplistic models (WLC) on the one hand, and accurate but expensive simulations (MC) on the other, still maintaining the analytical aspect which allows one to draw conclusions of physical interest. This is achieved using [41,42] as a starting point for bridging the two historical lines of research, i.e., exploiting efficient (real) path integral techniques in the semiclassical approximation [48][49][50][51][52][53] (or Laplace method [54]), and working within the special Cosserat theory of rods in SE (3). Namely, for studying the end-to-end relative displacements of a fluctuating polymer at thermodynamic equilibrium with a heat bath, we describe the configurations of the chain in a continuum limit by means of framed curves over the special euclidean group. Thus, from an assumed Boltzmann distribution on rod configurations, a conditional probability can be expressed as the ratio of a Boltzmann weighted integral over all paths satisfying the desired end conditions, to the analogous weighted integral over all admissible paths (partition function). The resulting path integrals are finally approximated via a quadratic (semiclassical) expansion about a minimal energy configuration, for which the crucial assumption is that the energy required to deform the system is large with respect to the temperature of the heat bath. This means computing probabilities for length scales of some persistence lengths or less, which turns out to be of great relevance in biology.
Although the present study is general and is applicable to various end-to-end statistics, we focus on the computation of ring-closure or cyclization probabilities for elastic rods, targeting three significant aspects. The first is the possibility of systematically distinguishing between the statistics provided by end positions alone (marginal looping) and the ones provided including also end orientations (full looping) [48], for Kirchhoff as well as for Cosserat rods. We emphasize that although Kirchhoff rod theory [55] generalises both Euler's elastica theory to model deformations in three dimensions, and the WLC model allowing arbitrary bending, twisting, and intrinsic shapes of the rod, it does not allow extension or shearing of the rod centerline. This is indeed a prerogative of the Cosserat, more general framework, where the centerline displacement and the cross-sectional rotation are considered as independent variables. We show that these additional degrees of freedom are crucial in the analysis of polymer chains in short-length scale, or stiff, regimes, in both the full and marginal cases, where the system exploits extension and shear deformations for minimizing the overall elastic energy, in the face of an increasingly penalizing bending contribution. This allows the cyclization probability density to take high values even when the WLC model (and Kirchhoff) is vanishing exponentially.
The second is addressing the "perfect problem" in the semiclassical context, where the symmetry of isotropy gives rise to a "Goldstone mode" [56] leading to a singular path integral, and requires a special treatment by suitably adapting (imaginary) quantum mechanical methods [57][58][59][60][61][62] and functional determinant theories [63][64][65], which are novel in such a generality in the context of elastic rod. For simple models, an analysis in this direction is present in [66]. The concepts of isotropy and nonisotropy can be roughly related to a circular shape rather than an elliptical shape for the cross section of the rod, and the two cases have two different mathematical descriptions in terms of Gaussian path integrals, which we discuss in detail in the course of this article. In particular, the effect of nonisotropy for semiflexible chain statistics has been addressed from a path integral point of view in [42] for the planar case and in [41] for the three-dimensional case (and will be here taken up and simplified), but without resolving the singularity arising in the isotropic limit.
The last significant aspect included in the present work is deriving approximated solution formulas that can always be easily evaluated through straightforward numerical solution of certain systems of Hamiltonian ODE, which in some particularly simple cases can even be evaluated completely explicitly. Versions of the solution formulas, involving evaluation of Jacobi fields at different equilibria and subject to different initial conditions (ICs), are obtained for the two cases of full and marginal ring-closure probabilities. The efficiency aspect in computing looping probabilities, maintaining the same accuracy of MC in the biologically important range less than one or two persistence lengths, is fundamental. This is because MC simulation is increasingly intractable due to the difficulty of obtaining sufficiently good sampling with decreasing polymer length, which is the limit where the approximation is increasingly accurate. Contrariwise our approximations are increasingly inaccurate in longer length regimes where good MC sampling is easily achieved. Remarkably, the qualitative behavior of the probability densities coming from Laplace approximation and from MC sampling are the same regardless of the length scale.
We stress that the stiffness parameters expressing the physical properties of the polymer are allowed to vary along the material parameter of the curve, leading to a nonuniform rod which, in the context of DNA, would represent sequencedependent variations. In addition, the model allows coupling between bend, twist, stretch and shear, as well as a nonstraight intrinsic shape. Notwithstanding the latter generality, we prefer to illustrate our method with some basic examples of uniform and intrinsically straight rods and comparing it with a suitable MC algorithm, in order to highlight the contributions provided by the different choices of cyclization boundary conditions (BCs) in the presence of isotropy or nonisotropy, and to investigate the effect of shear and extension when moving from Kirchhoff to Cosserat rods. Finally, the results will be FIG. 1. Schematic representation of a Cosserat rod with an elliptical cross section (nonisotropic), where bending, twist, shear, and extension/compression are allowed deformations. For Kirchhoff rods only bending and twist are permissible. The standard WLC model takes only bending deformations into account, and the rod is assumed to be intrinsically straight with circular cross section (isotropic). exposed under the hypothesis of linear elasticity, even though the theory applies to more general energy functionals.
The structure of the article is as follows. In Sec. II we give an overview of the statics of special Cosserat rods, with particular emphasis on equilibria and stability for the boundary value problems (BVPs) involved, and we further establish the relations with simpler models. In particular, the Hamiltonian formulation of the Euler-Lagrange and Jacobi equations provides a common theoretical framework for both Kirchhoff and Cosserat rods. In Sec. III we set out a preview of the examples that will be considered in the course of the article, namely, in the context of linear elasticity. Here we focus on the physical properties that characterize shearable and extensible-compressible polymers and explain how these degrees of freedom improve the understanding of the problem. Therefore, we study the minimizers of the energy, distinguishing between the nonisotropic and isotropic cases. The role of the continuous variational symmetries of isotropy and uniformity is explained. Before describing the computational setting in detail, we devote a section (Sec. IV) for summarizing the general formulas that we obtain for estimating end-to-end probabilities of fluctuating elastic rods as a proxy for interpreting the behavior of polymers in a thermal bath. Then we introduce the path integral formulation of the problem in Sec. V, prescribing an appropriate parametrization of the rotation group and giving the functional representations of full and marginal looping probability densities. Afterwards, the explicit approximated formulas for such densities are derived, initially in the case of isolated minimizers and thereafter in presence of nonisolation, for which a special theoretical analysis is performed. Moreover, in Sec. VI we provide a MC algorithm for stochastic elastic rods, exploited to benchmark our results. The examples are finally investigated from the point of view of cyclization probabilities in Sec. VII, with special focus on shear and extension contributions for Cosserat rods in the short-length scale regimes. Further discussion and conclusions follow.

II. BACKGROUND ON ELASTIC ROD EQUILIBRIA AND THEIR STABILITY
A comprehensive overview of the theory of elastic rods in the context of continuum mechanics can be found in [67]. In particular, we follow the specific notation and Hamiltonian formulations introduced in [68]. Briefly, a configuration of a Cosserat rod is a framed curve q(s) = (R(s), r(s)) ∈ SE (3) for each s ∈ [0, L], which may be bent, twisted, stretched, or sheared. The vector r(s) ∈ R 3 and the matrix R(s) ∈ SO(3) model, respectively, the rod centerline and the orientation of the material in the rod cross section via a triad of orthonormal directors {d i (s)} i=1,2,3 attached to the rod centerline, with respect to a fixed frame {e i } i=1,2,3 . As a matter of notation, the columns of the matrix R(s) in coordinates are given by the components of the vectors d j (s) in the fixed frame Fig. 1 we show a schematic representation of the the degrees of freedom allowed within the special Cosserat theory of rods in relation to other simpler models that will be outlined in the course of this section.
Strains are defined as u(s), v(s) where d i = u × d i , r = v, with u the Darboux vector and the prime denoting the derivative with respect to s. Sans-serif font is used to denote components in the director basis (e.g., u i = u · d i ), and we write u = (u 1 , u 2 , u 3 ), v = (v 1 , v 2 , v 3 ), etc. Physically, u 1 and u 2 represent the bending strains and u 3 the twist strain. Analogously, v 1 and v 2 are associated with transverse shearing, whereas v 3 with stretching or compression of the rod. In compact form, we have u × (s) = R(s) T R (s), v(s) = R(s) T r (s), where u × is the skew-symmetric matrix or cross product matrix of u having (1,2), (1,3), and (2,3) entries, respectively, equal to −u 3 , u 2 , and −u 1 .
The stresses m(s) and n(s) are defined as the resultant moment and force arising from averages of the stress field acting across the material cross section at r(s). In the absence of any distributed loading, at equilibrium the stresses satisfy the balance laws n = 0, m + r × n = 0. Equilibrium configurations can be found once constitutive relations are introduced, which we do in a way that facilitates the recovery of the inextensible, unshearable limit typically adopted in polymer physics.
Namely, we consider a pair of functions W, W * :R 3 × R 3 × [0, L] → R that (for each s ∈ [0, L]) are strictly convex, dual functions under Legendre transform in their first two arguments, and with 0 ∈ R 6 their unique global minimum. Ifû(s) andv(s) are the strains of the unique energy minimizing configurationq(s), then ∀ > 0, we introduce the Hamiltonian function H = W * (m, n; s) + m ·û + n ·v, and the constitutive relations are u = ∂H /∂m = W * 1 (m, n; s) +û, v = ∂H/∂n = W * 2 (m, n; s) +v, which can be inverted Note the use of the subscripts to denote partial derivatives with respect to the first or second argument. The standard case of linear constitutive relations arises when W * (x; s) = 1 2 x · R(s)x and W (y; s) = 1 2 y · P (s)y for x, y ∈ R 6 , where R 6×6 P −1 (s) = R(s) = R(s) T > 0, with P (s) a general nonuniform stiffness matrix and R(s) the corresponding compliance matrix. For each > 0 and given W , W * , we arrive at a well-defined Cosserat rod theory, where, e.g., the full potential energy of the system might include end-loading terms of the form λ · [r(L) − r(0)], λ ∈ R 3 .
The point of the above formulation is that the Hamiltonian and associated constitutive relations behave smoothly in the limit → 0, which imply the unshearability and inextensibility constraint on the strains v(s) =v(s), wherev(s) are prescribed. This is precisely a Kirchhoff rod model, abbreviated as (K), in contrast to (C) for Cosserat. However, the → 0 limit of the (C) Lagrangian is not smooth; rather the potential energy density for the (K) rod is the Legendre transform of W * (m, 0; s) + m ·û + n ·v w.r.t. m ∈ R 3 , or W (K) (u −û; s) − n ·v. In the case of linear elasticity for a (C) rod with P (s) = There is an extensive literature concerning the study of equilibria of a given elastic rod. Numerically this involves the solution of a two-point BVP, which can reasonably now be regarded as a straightforward well-understood procedure. Often coordinates on SO(3) are introduced and the resulting system of second-order Euler Lagrange equations associated with the potential energy is solved numerically. We adopt an Euler parameters (or quaternions) parametrization of SO(3), but solve the associated first-order canonical Hamiltonian system subject to appropriate (self-adjoint) two-point BCs, so that the inextensible, unshearable (K) rod is a simple smooth limit of the extensible, shearable (C) case.

A and K(s), B(s),
In this article we are primarily interested in the two specific BVPs, denoted respectively by (f) and (m): The BVP (f) arises in modeling looping in SE (3) including the particular case of cyclization where r L = 0 and R L = 1.
The BVP (m) arises in modeling looping in R 3 , where the value of R L is a variable left free, over which one marginalizes. In general, for rod two-point BVPs, equilibria with given BCs are nonunique. For isotropic or uniform rods, and for specific choices of r L and R L in (f) and (m), equilibria can arise in continuous isoenergetic families [69], a case of primary interest here. As we assume hyper-elastic constitutive relations with stability of rod equilibria can reasonably be discussed dependent on whether an equilibrium is a local minimum of the associated potential energy variational principle. For (C) rods classification of which equilibria are local minima has a standard and straightforward solution. The second variation δ 2 E is a quadratic functional of the perturbation field h = (δc, δt), where the sans-serif font q(s) = (c(s), t(s)) ∈ R 6 is a given parametrization of SE (3) for the configuration variable in the director basis, which will be specified later in the article, and reads as where P(s), C(s), and Q(s) are coefficient matrices in R 6×6 computed at any equilibrium. The Jacobi equations are the (second-order) system of Euler-Lagrange equations for Eq. (4), or equivalently the linearization of the original Euler-Lagrange equations for the potential energy variational principle. One then solves a 6 × 6 matrix valued system, namely, an initial value problem for the Jacobi equations with ICs coinciding with the ones given later in the article when computing probability densities from Jacobi fields (shooting towards s = 0, where in both (f) and (m) Dirichlet BCs are present; the case with Neumann BCs at both ends is more delicate [70]). Provided that the determinant of the matrix solution does not vanish in [0, L), then there is no conjugate point and the equilibrium is a local minimum [71][72][73].
As described fully in [74], the constrained case of (K) is more subtle and a theory dating back to Bolza for isoperimetrically constrained calculus of variations must be applied [75]. However, the Hamiltonian version of the Jacobi equations for rods (just like the Hamiltonian version of the Euler-Lagrange equilibrium equations) has a smooth limit as → 0, and the limit corresponds to the Hamiltonian formulation of the Bolza conjugate point conditions as described in [72]. The Jacobi equations in first-order Hamiltonian form are written as with the Hamiltonian skew-symmetric matrix J = 0 1 −1 0 ∈ R 12×12 , E(s) the symmetric matrix driving the system which will be detailed later, and M(s) ∈ R 6×6 the conjugate variable of the Jacobi fields H(s) under the Legendre transform.
In the following, we assume the existence and stability of the minimizers of the elastic energy (3) q f and q m satisfying the BCs (f) in Eq. (1) and (m) in Eq. (2). Note that the intrinsic configuration of the rodq is itself a minimizer (global) satisfying Stability of equilibria is not the focus of this article, but we will show that the volume of certain Jacobi fields, i.e., the actual (positive) value of a Jacobi determinant, plays a central role in the evaluation formula for the quadratic path integrals that arise in our Laplace approximations to looping probabilities.
The connection between Jacobi fields and quadratic imaginary path integrals is well known in the case that the coefficient matrix C(s) in the cross-terms in Eq. (4) vanishes (or is symmetric and so can be integrated away). By contrast, for elastic rods a nonsymmetric C(s) is typically present, and the approach of Papadopoulos [49] is required to evaluate the quadratic path integrals; as described in [41,42] a further Riccati transformation for the Papadopoulos solution formula is necessary to recover a Jacobi fields expression. Moreover, in [48] the latter studies are generalized for different choices of BCs on the paths, in particular for dealing with the partition function and solving the marginalized problem.
The main contributions of this article are to demonstrate that the approach of [41,42] for conditional probability densities can be extended in two ways. First, isolated equilibria to BVP (m) can be treated, in addition to the case of isolated equilibria to BVP (f), and second, the case of nonisolated equilibria of both BVPs (f) and (m) (as arises for isotropic rods) can be handled by appropriately generalizing a particular regularization procedure [64,65] within Forman's theorem in the field of functional determinants [63]. Furthermore, the underlying physical phenomena arising from the different cases are discussed and explained within some guiding examples. For a polymer, the questions we are trying to answer would be interpreted as follows: what is a good estimate of the probability of the end monomers coming into contact with each other? How is the latter value changing if we im-pose an orientation constraint on the binding site? How does the shape of the cross section (isotropic or nonisotropic) affect the statistics? And, finally, what happens if we deviate from the standard inextensible and unshearable model and incorporate shear and extension as possible deformations?

III. A PREVIEW OF THE EXAMPLES CONSIDERED
The method developed in the present article will be applied, as a fundamental example, to a linearly elastic, uniform, with diagonal stiffness matrix, intrinsically straight and untwisted rod [P (s) Neither intrinsic shear nor extension is present. Since we are primarily interested in ring-closure or cyclization probabilities, we look for minimizers of the energy satisfying the BCs reported in Eqs. (1) and (2) with r L = 0 and R L = 1.
First, we consider a nonisotropic rod (k 1 = k 2 ), further assuming w.o.l.o.g. that k 1 < k 2 . For the case of full looping (f), there exist two circular, untwisted, isolated minima q f lying on the y − z plane characterized by u f = (±2π/L, 0, 0) and v f = (0, 0, 1). In particular, the one having nonpositive y coordinate is given by r f (s) = L 2π (0, cos (2π s/L) − 1, sin (2π s/L)) and the rotation matrix R f (s) is a counterclockwise planar rotation about the x axis of an angle , n f = 0 and the energy is simply computed as E (q f ) = 2π 2 k 1 /L. We observe that these solutions are special for the fact of being the same both for (K) and (C) rods, which is not the case in general. By contrast, there are no simple analytical expressions for the two planar and untwisted teardrop shaped isolated minimizers q m involved in the marginal looping problem (m), and elliptic functions or numerics must be used. For example, in the (K) case, the rotation angle ϕ m (s) can be derived using elliptic functions in terms of the constant unknown force n m = (0, n2, n3) [55,76,77]. The qualitative shapes of the minimal energy configurations are reported in Fig present study their contributions will be neglected because of their higher elastic energy.
We continue the presentation with a brief stability analysis, showing that the circle and teardrop solutions are stable, with exceptions for the (C) rod in the limit of the undeformed length L going to zero, where bifurcations occur. For (C) rods, cyclization problems (f) and (m) always admit a "compressed" trivial solution q c , characterized by r c = 0, R c = 1, u c = 0, m c = 0, v c = 0, n c = (0, 0, −a 3 ) with energy E (q c ) = a 3 L/2, which starts to play an important role (this is not mentioned in [42]). In summary, for the full (C) case it exists L f > 0 such that the latter solution becomes stable and has lower energy than the circular minimizer q f if 0 < L < L f . In this regime the system will be mainly driven by the compressed solution (even if the circle remains stable). Moreover, for the marginal (C) case, it exists L m > 0 such that the stable teardrop solution q m ceases to exist in the interval 0 < L < L m , merging with the compressed solution which becomes stable. In both the cases, the above observations will have a strong impact on the trend of the estimated cyclization probability densities, that is confirmed by MC simulations. More precisely, analyzing the determinant of the associated Jacobi fields (5) [with ICs and matrix E(s) given later in Eq. (8) and Eqs. (B3)-(B5)] by means of conjugate point theory, we observe that the compressed solution is stable (i.e., a minimizer of the energy) in the range 0 < L < L f for the (f) case, and in 0 < L < L m for the (m) case, where L f = 2π/a 3 min{ √ k 1 a 2 , √ k 2 a 1 } and L m = L f /2. Moreover, as already mentioned, for full looping (f) there exist also circular solutions q f , which are stable for all L > 0, with energy 2π 2 k 1 /L. (This is true except for k 1 > k 3 , L < 2π √ (k 1 − k 3 )/a 1 , but in the present article we will not treat such an instability of the circular solution.) Note that if k1 k2 and a 1 = a 2 = a 3 , then For marginal looping (m), the teardrop solution q m is not present in the interval 0 < L < L m , transforming into the compressed solution which becomes stable. We show the bifurcation diagrams in Fig. 4 for a nonisotropic (C) rod. Observe that E (q m ) does not explode for small lengths, but instead reaches a maximum and decreases towards E (q c ). By contrast, for a (K) rod the circular and teardrop solutions exist and are stable for all L > 0, with energy diverging approaching L = 0, and no compressed solution is present.
In addition to the above statements, the isotropic case requires a more detailed analysis for the presence of a continuous symmetry. Namely, for a general linearly elastic (transversely) isotropic (C) rod defined by P (s) = diag{k 1 (s), k 2 (s), k 3 (s), a 1 (s), a 2 (s), a 3 (s)} with k 1 = k 2 , a 1 = a 2 andû 1 =û 2 =v 1 =v 2 = 0, it is known [69] that for cyclization BCs (f) in Eq. (1) and (m) in Eq. (2) the equilibria are nonisolated and form a manifold obtained, starting from a known solution, by a rigid rotation of the rod of an angle θ about the z axis and a subsequent rotation of the framing by an angle −θ about d 3 (s), for θ ∈ [0, 2π ) (register symmetry). As a consequence, in our particular examples, once selected, e.g., the nonisotropic solution lying in the y-z plane, y 0 and characterized by the configuration (R(s), r(s)), s ∈ [0, L], then we get an entire family of min- where Q θ is defined as the counterclockwise planar rotation matrix about the z axis of an angle θ ∈ [0, 2π ) (Fig. 2). As a side note for the (f) example, being the circular solutions the same for (K) and (C) rods, the isotropy symmetry arises even if a 1 = a 2 .
Furthermore, for a general linearly elastic uniform rod, for which the stiffness matrix P and the intrinsic strainŝ u,v are independent of s, another continuous symmetry is present for the cyclization BCs (f) in Eq. (1). In fact, starting from a known solution characterized by the configuration (R(s), r(s)), s ∈ [0, L], it is possible to obtain a family of equilibria parametrized by s * ∈ [0, L) in the following way: select s * ∈ [0, L), rigidly translate the rod by −r(s * ), reparametrize the rod using the parameter t ∈ [0, L] such that s = t + s * (mod L), rigidly rotate the rod about the origin In the present article we will deal with only one symmetry parameter, namely, θ ∈ [0, 2π ) associated with isotropic rods, where the presence of a family of minimizers translates into a zero mode ψ α (s; θ ) (α standing both for f and m) of the self-adjoint operator S α associated with the second variation (4), as will be discussed in due course. Therefore, the stability analysis reported in Fig. 4 is totally analogous for the isotropic case, except from the fact that an entire family of minimizers is involved and a conjugate point is always present due to the zero mode. Furthermore, the theory can be applied to the uniformity symmetry alone and generalized to cases in which isotropy and uniformity allow the coexistence of two nondegenerate symmetry parameters (θ, s * ) generating a manifold of equilibria isomorphic to a torus, as it is the case of figure-eight minimizers with (f) cyclization BCs. Finally, note that in the following theory there is no assumption either of uniformity of the rod, or, in general, of a straight intrinsic shape.

IV. STATEMENT OF THE PROBLEM AND GENERAL RESULTS
In this section we describe the problem at the heart of this paper and present the general formulas that we derive in the context of end-to-end probabilities for fluctuating elastic rods, valid in both the (C) and (K) cases. The proof and the application of these results will follow in separate sections. Thus we consider an elastic rod at thermodynamic equilibrium with a heat bath in absence of external forces, assuming w.o.l.o.g. that q(0) = q 0 = (1, 0). Then, given a prescribed q L = (R L , r L ) ∈ SE (3), we formulate the problem of computing a conditional probability density function (pdf) for the other end of the rod to satisfy at s = L either q(L) = q L , or the weaker condition r(L) = r L . The first case gives rise to a conditional pdf (f) over the space SE (3) denoted by ρ f (q L , L|q 0 , 0), whereas the second one represents the R 3valued marginal (m) over the final rotation variable, with no displacement constraint on R(L), that will be denoted by ρ m (r L , L|q 0 , 0). The following results are given for the case of linear elasticity, although the theory developed in the article is general.
We show that an approximate form of the conditional probability density in the case of an isolated minimizer q α (s) of the elastic energy (3) [with respect to the associated BVPs (f) and (m)] reads as with We further show that an approximate form of the conditional probability density in the case of nonisolated minimizers q α (s; θ ), obtained by means of a suitable regularization procedure, reads as and we are interested in the cyclization values ρ f (q 0 , L|q 0 , 0), ρ m (0, L|q 0 , 0). In particular, μ ψ α ∈ R 6 and H α ∈ R 6×6 are, respectively, the conjugate momentum of the zero mode and the Jacobi fields associated with S α , both computed by means of Eq. (5) where χ is an arbitrary matrix with unit determinant such that the ith column corresponds to μ ψ f (L) and X = X 1,1 X 1,2 X 2,1 X 2,2 ∈ R 6×6 , partitioned in 3 × 3 blocks, is an arbitrary matrix with determinant equal to −1 such that the ith column corresponds to ([ψ m ] 1:3 , [μ ψ m ] 4:6 ) T (L).

V. FLUCTUATING ELASTIC RODS AND THE PATH INTEGRAL FORMULATION
If a polymer interacts with a solvent heat bath, the induced thermal motion gives rise to a stochastic equilibrium that we model making use of a Boltzmann distribution on rod configurations satisfying q(0) = q 0 [41,42], of the form Z −1 e −βE (q(s)) , with β the inverse temperature and Z the partition function of the system. A precise treatment of the previous expression requires the introduction of the path integral formalism [21,[50][51][52]. Namely, the SE (3) and R 3 densities ρ f and ρ m are respectively given as the ratios of infinite dimensional Wiener integrals [48]: The limits of integration are dictated by the BCs (1) and (2), respectively, and Z is a path integral over all paths with BCs given in Eq. (6) that guarantees the normalization condition: The prescriptions m(L) = 0 for K m and m(L) = n(L) = 0 for Z account for Neumann natural BCs at s = L and concern the minimizers. We stress that it is key that at this stage the model is an extensible, shearable rod, namely, with (C) energy (3), otherwise the problem could not be expressed as simple BCs at s = 0 and s = L. Moreover, to apply all the path integral machinery, we first have to deal with the rotation group SO(3), being part of the configuration variable q(s) = (R(s), r(s)), which gives rise to a manifold structure that should be treated carefully in order to recover eventually a "flat space" formulation.
Following [41], we show in Appendix A how to build an R 6 parametrization of SE (3) [SO (3) is not simply connected and a zero measure set of rotations is neglected] adapted to a given unit quaternionγ ∈ R 4 . In particular we make use of the Haar measure on SO(3) and derive the metric tensor associated with the parametrization. Namely, theγ-adapted parametrization of SE (3) denoted by q(s) = (c(s), t(s)) ∈ R 6 exploits the relation between unit quaternions (or Euler parameters) γ and elements of SO(3) and is given by (14) with c = (c 1 , c 2 , c 3 ) ∈ R 3 , R(γ ) the rotation matrix expressed byγ, and B 1 , B 2 , B 3 in R 4×4 reported in Eq. (A1). Moreover, by means of the Feynman discrete interpretation of the path integral measure [21], the metric tensor and the infinitesimal volume measure read, respectively, The latter results are implemented by choosing three different curves of unit quaternionsγ (s) to be the curves defined by the rotation component R(γ ) of the minimizers q f , q m andq respectively, which characterize the three different parametrizations involved in the computation of K f , K m , and Z in view of the semiclassical approximation. Then, replacing the configuration variable q(s) ∈ SE (3) with the sans-serif fonts q(s) ∈ R 6 , we can formally write the integrand and measure in Eqs. (12) and (13) as e −βE (q) √ det [g(c)] Dq. The treatment of the metric factor relies on the introduction of real-valued ghost fields for exponentiating the measure, as can be found in [58]. This means rewriting the factor as a Gaussian path integral in the ghost field z(s) ∈ R 3 satisfying z(0) = 0 with energy 1 2 L 0 z T g −1 (c)z ds. After that, we consider the path integral expressions in the joint variable w = (q, z), e.g., In the following, even if the theory could be given in principle for a general strain energy density W , in order to perform concrete computations we refer to the case of linear elasticity, where W is a quadratic function, driven by the stiffness matrix P (s): Moreover, we also refer to the particular looping case of ring-closure or cyclization, evaluating ρ f at q L = q 0 and the marginal ρ m at r L = 0; the same conditions apply to the minimizers.

A. Looping probabilities in the case of isolated minimizers
Since the elastic energy functional (17) is nonquadratic in q, after the parametrization we approximate K f , K m , and Z by means of a second-order expansion about a minimal energy configuration [48][49][50][51][52][53], known as the semiclassical method, or, in our real-valued context, Laplace expansion [54]. The present work follows the setup of [48]. We further recall that such an approximation holds when the energy required to deform the system is large with respect to the temperature of the heat bath, i.e., in the short-length scale, or stiff, regimes.
First, note that there is no contribution to the result coming from the ghost energy when approximating path integrals of the kind of Eq. (16) to second order in the joint variable w. This is a consequence of the structure of the metric tensor (15), i.e., g −1 (c) = (1 + c · c)(1 + c ⊗ c), and therefore we can consider only the elastic energy (17)  In the present case of linear elasticity, the second variation (4) is characterized by P, related to the stiffness matrix P, and C, Q which can be computed as follows in terms of strains, forces, and moments of the minimizer involved, generically denoted byq = (R(γ ),r). In elastic rod theory, the natural parametrization for the variation field aroundq is directly provided by the Lie algebra so(3) of the rotation group in the director frame, namely, δR = R(γ )δη × , where δη × denotes the skew-symmetric matrix or cross-product matrix of δη ∈ R 3 . In order to show the relation between δη and the variation field δc, we use the formula δη = 2( 3 i=1 e i ⊗ B iγ )δγ (which is substantially the relation between the Darboux vector and Euler parameters; see, e.g., [68]) with δγ = ∂γ ∂c | c=0 δc = 3 j=1 B jγ δc j referring to Eqs. (A1) and (14), and we conclude that δη(s) = 2 δc(s).
With reference to [78], the second variation of the linear hyper-elastic energy (17) in the director variable ω = (δη, δt) where P is the stiffness matrix and C, Q are respectively given in terms of strains, forces, and moments by Eqs. (B1) and (B2). Finally, introducing the matrix D = 21 0 0 1 , we have that the second variation in the variable h = (δc, δt) (4) is given by The Jacobi equations in first-order Hamiltonian form associated with the latter second variation functional are given in Eq. (5) and are driven by the symmetric matrix E(s) ∈ R 12×12 detailed in Eq. (B3). The Jacobi fields H(s) ∈ R 6×6 , together with the conjugate variable under the Legendre transform M(s) ∈ R 6×6 represent the solutions of the Jacobi equations once prescribed appropriate ICs. The columns h of H and the ones μ of M are related by μ = Ph + Ch.
Note that until now the formulation adopted is for the general (C) rod with extension, shear, and hence an invertible stiffness matrix P. The constrained inextensible and unshearable case (K) requires the stiffness components B and A to diverge (as discussed in [41,42,68]), specifically as B/ and A/ 2 , for → 0. Switching to the Hamiltonian formulation, given a (C) rod the compliance matrix R (which is the inverse of P) has a smooth limit for → 0. Namely, for a (K) rod we recover R(s) = R 1,1 In conclusion, once prescribed a symmetric and positive definite matrix K (K) = R −1 1,1 , there exists a sequence of positive definite and symmetric compliance matrices for the (C) case converging smoothly to the (K) case, implying that the expressions (B4) and (B5) for the blocks of the matrix E(s) (B3) hold for both (C) and (K) rods. We emphasize that for the (K) case δ ∂W ∂t is a basic unknown of the Jacobi equations and cannot be found using the relation μ = Ph + Ch, since the latter is not defined.
The resulting path integrals arising from the semiclassical method are of the form, e.g., and similarly for K m and Z but considering the different minimizers and linearized BCs. Then, applying the results derived in [48] for Gaussian path integrals, which are in turn extensions of [49], we recover the approximate form of the conditional probability density (7). In principle, denoting bŷ H(s) ∈ R 6×6 the Jacobi fields computed atq subject to the ICsĤ(L) = 1,M(L) = 0 [48], the numerator and denominator in Eq. (7) should be, respectively, e −β(E (q α )−E (q)) and det [H αĤ −1 (0)], in order to include the contribution coming from the evaluation of the partition function Z. However, the result simplifies since E (q) = 0, beingq the intrinsic configuration of the rod. At the same time E 1,1 is the zero matrix for this case, which impliesM(s) = 0 ∀s [according to the IĈ M(L) = 0] and consequentlyĤ(s) must satisfy a linear system whose matrix has zero trace. Thus, by application of the generalized Abel's identity or Liouville's formula, ∀s we have that det[Ĥ(s)] = det[Ĥ(L)] = det[1] = 1. Furthermore, it is worth to mention that here the partition function computation is not affected by approximations, even if it apparently undergoes the semiclassical expansion. In fact, there exists a change of variables presented in [41,42] which allows an equivalent exact computation exploiting the specific BCs involved in Z. In general, the latter change of variables is not applicable and the present method must be used, e.g., for nonlinear elasticity or in the case of a linearly elastic polymer subject to external end loadings, for which the shape of the energy leads to a nontrivial contribution of the partition function that must be approximated.

B. Looping probabilities in the case of nonisolated minimizers
In this section we consider nonisolated minimizers arising as a consequence of continuous symmetries of the problem. In particular, we provide a theory for one symmetry parameter, namely θ ∈ [0, 2π ) (as we want to deal with isotropic rods), but the same scheme can be suitably generalized to more symmetry parameters. The presence of a family of minimizers denoted by q α (s; θ ) translates into a zero mode ψ α (s; θ ) = ∂ ∂θ q α (s; θ ) [53] of the self-adjoint operator S = −P d 2 ds 2 + (C T − C − P ) d ds + Q − C associated with the second variation (4), namely, δ 2 E = (h, Sh), where (·, ·) is the scalar product in the space of square-integrable functions L 2 ([0, L]; R 6 ). Consequently, we cannot proceed as before, for otherwise expression (7) will diverge for the existence of a conjugate point at s = 0.
Thus, in evaluating expression (12) for K f and K m , we adapt the parametrization to the minimizer corresponding to θ = 0, our choice of the gauge in applying the collective coordinates method, which amounts to a Faddeev-Popov-type procedure [57], widely used in the context of quantum mechanics for solitons or instantons [59][60][61][62], of inserting the Dirac δ transformation identity within the path integral, in order to integrate over variations which are orthogonal to the zero mode. Once performed the semiclassical expansion as before about q α , exchanged the order of integration Dh ↔ dθ to get a contribution of 2π , and having approximated to leading order both the metric tensor and the factor |∂/∂θF | θ=0 ≈ ψ α (s; 0) , we are left with the computation of a ratio of Gaussian path integrals for the linearized parametrized BCs associated with Eq. (1) and Eq. (2), respectively. For notation simplicity, throughout this section S stands for β 2π S andŜ for β 2πŜ , the latter operator driving the Gaussian path integral Z g arising from the partition function Z, in which the minimizerq is isolated. Note that, since the argument of the delta distribution must vanish for θ = 0 according to Eq. (19), then the integration for the numerator is performed on the minimizer q α (s; 0) with associated zero mode ψ α (s; 0); in the following they will both be denoted simply by q α and ψ α .
Interpreting Eq. (20) as Det(Ŝ)/Det (S), i.e., the square root of the ratio of the functional determinants for the oper-atorsŜ and S, the latter with removed zero eigenvalue (thus the star symbol) [64,65], we consider the following general strategy for its evaluation. Given the second variation operator S acting on h(s) ∈ R 6 , with s ∈ [0, L] and BCs determined by the square matrices T 0 and T L as T 0 μ(L) = 0, we state Forman's theorem [63] in Hamiltonian form as for W(s) ∈ R 12×12 whose columns (h, μ) T solve the homogeneous problem Sh = 0 (i.e., the Jacobi equations (5) with the extra β 2π factor, completed as W = JEW), and the trivial partition function contribution has already been evaluated. It is important to note the freedom of choosing W(0), W(L) consistently; the latter statements are justified by the following considerations.
Given two matrix differential operators Ω = G 0 (s) d 2 ∈ [a, b], the results of Forman [63] provide a simple way of computing the ratio of functional determinants Det(Ω)/Det(Ω), once prescribed the BCs I a The idea is now to compute expression (21) for the operator S subject to carefully chosen perturbed BCs T (ε) 0 , in order to avoid the zero mode. This gives rise to a quasizero eigenvalue that can be found analytically using our extension to general second variation operators (including cross-terms) of the trick introduced in [64]. Finally, by taking the limit for ε → 0 in the ratio of the regularized expression (21) to the regularized quasizero eigenvalue, we recover the desired quantity Det (S)/Det(Ŝ).
We anticipate here the results for the approximation formulas of the probability densities in the case of nonisolated minimizers (already stated in Eqs. (9) and (10) when presenting the final formulas), valid also for (K) rods as detailed in the previous section (note that the factor ψ α simplifies out within the regularization procedure) and we are interested in the cyclization values ρ f (q 0 , L|q 0 , 0), ρ m (0, L|q 0 , 0). In particular, μ ψ α ∈ R 6 and H α ∈ R 6×6 are, respectively, the conjugate momentum of the zero mode and the Jacobi fields associated with S α , both computed by means of Eq. (5) where χ is an arbitrary matrix with unit determinant such that the ith column corresponds to μ ψ f (L) and X = X 1,1 X 1,2 X 2,1 X 2,2 ∈ R 6×6 , partitioned in 3 × 3 blocks, is an arbitrary matrix with determinant equal to −1 such that the ith column corresponds to ([ψ m ] 1:3 , [μ ψ m ] 4:6 ) T (L).
We are now ready to explain how to regularize the functional determinants for S f and S m respectively, in order to get rid of the zero eigenvalue. Starting from the pure Dirichlet case, the BCs are given as T The last step consists of finding the nonzero eigenvalue λ α (ε) associated with the eigenfunction ψ α (ε) (arising from the zero mode ψ α ) of the operator S α with perturbed BCs. First, we have that (ψ α , S α ψ α (ε) ) = λ α (ε) (ψ α , ψ α (ε) ), and the left-hand side can be rewritten as (ψ α , where the second equality comes after integration by parts and the third and fourth ones are a consequence of the BCs. Finally, being λ α We conclude with a technical remark. We observe that a priori the solution formulas for isolated minimizers could be recovered by applying Forman's theorem in the framework of functional determinants (as done here for the nonisolated case); however, there we exploit the insightful connection with the more standard theory of path integrals via "time slicing." Exploring the latter possibility not only allows us to gain a deeper understanding of the subject, but is crucial to developing the right ideas for using Forman's formalism in a more general setting.

VI. A MONTE CARLO ALGORITHM FOR STOCHASTIC ELASTIC RODS
In this section we refer to the approach of [45][46][47] for DNA MC simulations of J-factors, using the "half-molecule" technique [44] for enhancing the efficiency. Namely, we give a Monte Carlo sampling algorithm for fluctuating linearly elastic rods according to the Boltzmann distribution having partition function (13), i.e., Z = q(0)=q 0 e −βE (q) Dq with energy (17), and we use the compact notation u = u −û, v = v −v for the shifted strains. First of all, we need to rewrite the infinite-dimensional problem as a finitedimensional one by means of a "parameter slicing method." This is achieved, after parametrizing the configuration variable as q(s) = (c(s), t(s)) ∈ R 6 , setting = L n with n a large positive integer and s j = j for j = 0, . . . , n. Moreover, by exploiting the change of variables (c j , t j ) → (u j , v j ) as presented in [41], we get the following equality up to a constant factor for the discrete version of the partition function Z: with 2 and the subscript j indicates that the associated term is evaluated in s j . We observe that the Jacobian factor J can be neglected, as discussed in [79], leading to the Gaussian distribution j=0 du j dv j which can be easily sampled by a direct MC method in order to get random instances of u j , v j , j = 0, . . . , n − 1, associated with a random framed curve with initial data q 0 = (1, 0). Note that, in the proposed uniform example with diagonal stiffness matrix, the Gaussian factorizes and the sampling is simply performed componentwise in terms of independent univariate Gaussians.
Since the conditional probability density is a function of the variables R L , r L , we need to reconstruct R n , r n from the sampled strains by discretization of the differential equations γ (s) = 1 R(γ (s))v(s), with [u(s)] i the ith component of u and R(γ ) the rotation matrix associated with the quaternion γ. This is achieved, e.g., by application of the scalar factor method, derived in [80] and discussed in [81], which is an efficient and precise one-step method for integrating the Darboux vector u, preserving the unit norm of the quaternion. Defining subject to the initial data γ 0 = (0, 0, 0, 1), and consequently r j+1 = r j + R(γ j )v j , r 0 = 0.
In the spirit of [47] for computing cyclization densities, we are now able to generate MC trajectories and assess whether or not q n = (R n = R(γ n ), r n ) is falling inside the given small region R ζ ,ξ of SE (3) centered in (1, 0) parametrized as the Cartesian product B ζ × B ξ of two open balls in R 3 , centered in 0, of radius ζ , ξ > 0 respectively. Namely, (R n , r n ) ∈ R ζ ,ξ if and only if c(γ n ) < ζ and r n < ξ, with c ∈ R 3 the same parametrization of SO(3) presented above, adapted toγ = (0, 0, 0, 1). Note that, since c(γ n ) = [γ n ] −1 4 ([γ n ] 1 , [γ n ] 2 , [γ n ] 3 ) and γ n = 1, the condition c(γ n ) < ζ is equivalent to [γ n ] −2 4 − 1 < ζ . Moreover, we have the following link between the probability of the set R ζ ,ξ [P(R ζ ,ξ )] computed using MC simulations and the conditional probability density defined in the theoretical framework where the notation | · | stands for the number of elements of a discrete set or the measure of a continuous set, and the accuracy of the approximation increases with n → ∞, |{all samples}| → ∞, ζ → 0, ξ → 0. The set R ζ ,ξ is measured by means of the product of the Haar measure and the Lebesgue measure for the SO(3) and the E (3) components. Thus, making use of the parametrization, Regarding the marginal ρ m (r 0 , L|q 0 , 0), the method is applied only considering the condition on r n for being inside the open ball B ξ with measure |B ξ | = 4πξ 3 /3, and neglecting all the details concerning the rotation component. More specifically, in order to enhance the efficiency of the algorithm, we refer to the approach adopted in [45][46][47] for DNA MC simulations, using the "half-molecule" technique as developed by Alexandrowicz [44]. In this technique, one computes M random instances each of the first and second halves of the framed curve and then considers all first-half-secondhalf pairs in order to generate M 2 random curves, allowing a large sample size contributing for each density data point and providing the necessary accuracy to the estimation. In particular, we give here the specifications for the simulations reported in the following section. For the (f) computations, ∼10 15 samples were produced for each data point, choosing n = 200 and ζ , ξ ranging from 2.5% to 6.6% of the parameter L. The estimated density value corresponds to the mean taken over 81 "boxes", along with the standard deviation for these boxes defining the range of the bar for each MC data point. For the (m) cases, ∼10 13 samples were produced for each data point, choosing n = 200 and ξ ranging from 0.1% to 4% of the parameter L; 40 different "boxes" were used for the final estimation.

VII. RESULTS AND DISCUSSION FOR THE EXAMPLES CONSIDERED
This section is dedicated to the application of formulas (7) and (9) in order to predict cyclization probabilities in a concrete example of a fluctuating polymer modeled as a linearly elastic, uniform, with diagonal stiffness matrix, intrinsically straight and untwisted rod [P (s) = P = diag{k 1 , k 2 , k 3 , a 1 , a 2 , a 3 },û = 0,v = (0, 0, 1)], as presented above. The chosen example allows the physical peculiarities of the problem to be investigated in a clear and effective manner, while also providing analytical expressions for particularly simple cases and capturing the phenomena involved. We remark that the theory proposed in this article is general and can be applied to nonuniform problems, e.g., to consider sequence-dependent variations in stiffness in the context of DNA modeling, as well as sequence-dependent intrinsic curvature.
We start with a preliminary analysis. Since in the (C) case the compressed (isolated) solution is a minimizer for the shortlength scale regimes, we evaluate analytically its contribution ρ c α to the cyclization probability density (f) and (m) for 0 < L < L f and 0 < L < L m respectively. Making use of Eq. (7) with ICs (8) and setting the nondimensional lengthL = L/l p for a given l p > 0, we get where x = x(α) with x( f ) = 1, x(m) = 2 and E p = β l p a 3 /2, ϑ 1 = (l p a 3 )/(2 √ k 1 a 2 ), ϑ 2 = (l p a 3 )/(2 √ k 2 a 1 ), with The latter formula is valid both for isotropic (setting k 1 = k 2 , a 1 = a 2 ) and nonisotropic rods. In the following we focus on the contribution ρ α to the cyclization probability density (f) and (m) coming from the circular and teardrop minimizers, respectively.

A. Nonisotropic polymers
First, we consider a nonisotropic rod (k 1 = k 2 ), further assuming w.o.l.o.g. that k 1 < k 2 . For the case of full looping (f), there exist two circular, untwisted, isolated minimizers q f lying on the y-z plane with energy 2π 2 k 1 /L. The existence of a couple of reflected minima simply translates into a factor of 2 in front of Eq. (7) and the semiclassical expansion is performed about one of them (e.g., about the one having nonpositive y coordinate). For this case Eq. (5) is a constant coefficients Jacobi system, that we solve analytically together with the first set of ICs in Eq. (8), in order to obtain the approximated formula for the cyclization probability density ρ f (q 0 , L|q 0 , 0) for both (C) and (K) rods. Setting the length scale l p = 2βk 1 , which corresponds to the planar tangenttangent persistence length for the same rod but constrained in two dimensions [79], and the nondimensional lengthL = L/l p , we get ρ f ≈ 2 e −π 2 /L h I h O , where h I and h O are the in-plane (of the minimizer) and out-of-plane contributions , ω 1 = k 3 /(a 1 l 2 p ). The (K) case is recovered setting a = b = 1, and the density obtained disregarding the factor h O coincides with the cyclization probability density for planar rods given in [42]. Note that the in-plane and out-of-plane contributions are computed by performing two separated Gaussian path integrals for the in-plane and out-of-plane variation fields, exploiting the decomposition of the second variation in two distinguished terms [41]. Moreover, the expressions in Eqs. (28) and (29) are valid under stability assumptions for k 1 < k 2 , k 1 = k 3 and equal to the limit k 3 → k 1 , i.e., ν 3 → 1, if k 3 = k 1 . We further underline that Eq. (29) diverges in the isotropic limit k 2 → k 1 , i.e., ν 2 → 1. The results for the full looping conditional probability in the case of isolated minimizers given above were first derived in [41], where the Gaussian path integrals are carried out in the variables ω = (δη, δt) instead of h = (δc, δt) as done here. As a consequence of the latter choice, in [41] all the formulas have a factor of 8 in front corresponding to the Jacobian factor of the transformation (actually in the cited work a factor of 2 is present, but it is a typographical error, it should be 8). We cite this reference for the explicit evaluation of the Jacobi fields leading to Eqs. (28) and (29).
In general, for computing the density ρ m (0, L|q 0 , 0) from Eq. (7) together with the second set of ICs in Eq. (8), numerics must be used. In fact, for the case of marginal looping (m), there are no simple analytical expressions for the two planar (y-z plane) and untwisted teardrop shaped isolated minimizers q m . However, in the (K) case there exists a scaling argument in the variable L, which allows one to provide a qualitative expression. Namely, given the fact that we can compute numerically a (K) equilibrium q m p for a given rod length l p , characterized by r p (s p ), where E p and h p have to be computed numerically, and the factor 2 accounts for the contribution of both the minimizers. By contrast, a simple scaling argument is not present for a (C) rod, therefore allowing for more complex behaviors. We show the results in Fig. 5 for a specific choice of the parameters, in the range L > L f and L > L m , respectively, for (f) and (m), so that the only accounted minimizers for the computation of the cyclization probability densities are the circular and the teardrop solutions, and we can apply Eqs. (29) and (30). The simulations show good agreement between the Laplace approximation and MC in the target small length domain. Even though the second-order expansion loses its quantitative power for larger lengths, the qualitative behavior is captured and the error does not explode. We recall that looping is a rare event and MC simulations are usually expensive and unfeasible; by contrast, the method proposed in the present article is performing successfully with much higher efficiency. It is also important to underline that for the specific example considered the difference in ρ f between (K) and (C) rods is due only to Jacobi fields, since the energy factor is the same, the circular minima having no extension and no shear deformations. The marginal case (m) is more representative of the general behavior where (K) and (C) minimizers are distinct solutions, which is true also for (f) BCs for arbitrary (nonuniform, with nonstraight intrinsic shape) elastic rods. In fact, in the short-length scale regimes, the possibility to exploit the additional degrees of freedom associated with extension and shear is crucial for minimizing the overall elastic energy, in the face of an increasingly penalizing bending contribution. This phenomenon allows the probability density to be remarkably higher than the (K) case below the persistence length, remaining almost constant and even increasing in the range where for the (K) rod (and therefore also for the WLC model) is exponentially vanishing. By contrast, for large lengths extension and shear become negligible. In addition, as a general statement, the Jacobi factor is fundamental to determine the peak of the density, in a domain where the energy is monotonically decreasing with length. On the other hand, the energy contri-FIG. 5. Comparison of cyclization densities between the path integral (PI) Laplace approximation (continuous lines) and MC (discrete points with standard deviation error bars) for a nonisotropic rod. For (K) we set β = 1, k 1 = 0.5, k 2 = 5, k 3 = 10; for the (C) case, we also set a 1 = a 2 = a 3 = 100. The quantities are reported in nondimensional form. In particular, the undeformed length of the rod is expressed in units of real persistence length l p ≈ 0.9, the harmonic average of k 1 and k 2 . In panel (a) we address the (f) case, reporting the values for ρ f and displaying in red the zero-order contribution. In panel (b) the results for the marginal density ρ m are reported, with a zoom window in log 10 scale in order to underline the peculiar small length trend. In this case, (K) and (C) rods differ in zero-order contribution of the energy, and two different red curves are displayed. bution dominates the system for smaller lengths. Finally, we clearly observe overall higher values for the marginal density compared to the full case because of the less restrictive BCs.

B. Isotropic polymers
Now we consider the isotropic case, i.e., k 1 = k 2 , a 1 = a 2 and a one-parameter family of nonisolated circular or teardrop minima arises. Given the minimizer in the y-z plane with y 0 represented by r(s) = (0, r 2 (s), r 3 (s)) and R(s) a counterclockwise planar rotation about the x axis of an angle ϕ(s), the one-parameter family of minimizers can be expressed as r(s; θ ) = Q θ r(s) and R(s; θ ) = Q θ R(s)Q T θ , where Q θ is defined as the counterclockwise planar rotation about the z axis of an angle θ ∈ [0, 2π ). Thus, taking the derivative of such minimizers with respect to θ and finally setting θ = 0, the zero mode can be easily recovered in the chosen parametrization to be ψ(s) = (0, 1 2 sin (ϕ), 1 2 [cos (ϕ) − 1], −r 2 , 0, 0). Moreover, the conjugate momentum of ψ is derived in general (for both (C) and (K) rods) substituting the zero mode itself and its to-be-found moment as the unknowns of the Jacobi equations in Hamiltonian form (5) computed on the minimizer associated with θ = 0 (recalling to multiply E 1,1 by β 2π and E 2,2 by 2π β ), and reads as μ ψ (s) = (0, βk 1 2π [cos (ϕ) + 1]ϕ , − βk 1 2π sin (ϕ)ϕ , − β 2π n 2 , 0, 0). At this point it is straightforward to apply the theory developed above for nonisolated minimizers, choosing χ to be a matrix with unit determinant such that the second column (i = 2) corresponds to μ ψ f (L) and X a matrix with determinant equal to −1 such that the fourth column (i = 4) corresponds to ([ψ m ] 1:3 , [μ ψ m ] 4:6 ) T (L), according to Eq. (10). Consequently, the ICs for the Jacobi equations are well defined, the energy is computed, e.g., for the minimizer corresponding to θ = 0 as before, and Eq. (9) is analytical for with a = 1 + (2π/L) 2 (η 1 + η 3 ) and all the other quantities have been defined previously. In particular, h I is the same as for the nonisotropic case and therefore the zero mode arises for the out-of-plane factor for which the above regularization is applied. The (K) limit is recovered as before setting For the marginal density ρ m numerics must be used, but in the (K) case we can carry on the scaling argument in the variable L as before, obtaining for given It is interesting to note that formulas (29)-(32) scale differently with length as far as the second-order correction term is concerned. The latter scalings naturally arise from the ones observed within simpler WLC models; see chapter 7 in [23]. The comparison between Laplace and MC simulations for isotropic polymers is shown in Figs. 6(a) and 6(b), for the same parameters addressed in the nonisotropic case, but now sending k 2 → k 1 . Once more time we only consider the FIG. 6. Comparison of cyclization densities between the path integral (PI) approximation and MC for an isotropic rod. For (K) we set β = 1, k 1 = k 2 = 0.5, k 3 = 10; for the (C) case, we also set a 1 = a 2 = a 3 = 100. The quantities are reported in nondimensional form, and the undeformed length of the rod is expressed in units of real persistence length l p ≈ 0.5. In panels (a) and (c) we address the (f) case, reporting the values for ρ f and displaying in red the zero-order contribution. The behavior for (C) in the small-length regime is shown in (c). In panels (b) and (d) the results for the marginal density ρ m are reported, with a zoom window in log 10 scale; the two different zero-order contributions for (K) and (C) rods are displayed in red. The behavior for (C) in the small-length regime is shown in panel (d).
contributions of the manifolds made of circular and teardrop minimizers, setting L > L f and L > L m . The fact that now k 2 is ten times smaller than the same parameter adopted in Fig. 5 implies that the overall trend of the density is shifted to the right in units of persistence length, allowing large effects of shear and extension compared to the more standard inextensible and unshearable models, as already discussed. We further observe that the approximation error is generally higher for (C) rods and for marginal looping (m), which is a consequence of the semiclassical expansion that depends on the stiffness values and BCs. For the simple examples considered, there clearly exist more accurate formulas for the (K) case in the literature; e.g., Eq. (32) can be related to the WLC formula (7.68) (p. 266 in [23]). However, the power of the method explained above lies in its generality and ability to easily provide approximation formulas for a wide range of potentially realistic and complex problems in the short-length scale regimes. By contrast, since the (C) case represents itself a novelty, we believe that basic examples are still important to understand the underlying physical behavior.
It is natural to ask what happens for L L f and L L m , respectively, in (f) and (m), for (C) rods [for (K) the former analysis based only on circular and teardrop solutions is valid for all lengths]. Due to the presence of the stable compressed solution in this range, the density diverges for vanishing length, and this is true for both isotropic and nonisotropic rods. In particular, for (f) here we sum up the contributions coming from the compressed solution (27) and the manifold of circular minimizers (31); for (m) only the compressed solution is present and we apply Eq. (27). At the critical lengths L f and L m a conjugate point arises for the compressed solution (in (m) the conjugate point arises also in the teardrop minimizer) and the Jacobi fields are singular, leading to an incorrect explosion of the probability density, which should be regularized. We do not address such regularizations, but in Figs. 6(c) and 6(d), we report the results for this length regime, together with MC simulations which connect our approximation formulas valid on the left and on the right of the singularities. We remark that the diverging behavior of the conditional probability density at zero length observed for Cosserat polymers is a consequence of the linearly elastic hypothesis on the energy functional and cannot be regarded as a physical behavior in the context of polymers made of discrete elementary units. However, we show the existence of a length range, not affected by a compressed stable solution, where high looping probabilities occur due to an energy relaxation of the minimizers achieved by exploiting the degrees of freedom associated with extension and shear.
Finally, in order to highlight the effect of shear and extension for larger lengths, in Fig. 7 we compare the (K) and the (C) cases in terms of the length and the value of the probability density at which the maximum of ρ f occurs, the first increasing and the second decreasing in presence of extension and shear.

VIII. CONCLUSIONS
In the present article we addressed the problem of computing looping probabilities from a continuum perspective, for different choices of BCs, with particular emphasis on extensible and shearable polymers, which are not generally treated in the standard literature of WLC-type models. Moreover, the proposed theoretical framework employed for deriving general looping formulas is supplemented with concrete examples, the results of which are also supported by extensive Monte Carlo simulations.
In a first approximation DNA fits the WLC hypothesis of inextensibility and unshearability. However, contradictory results have been reported for DNA below the persistence length since the studies of Cloutier and Widom [82], actually showing enhanced cyclization of short DNA molecules not explainable by WLC-type models. In a recent study [31] the authors conclude that "determining whether the high bendability of DNA at short-length scales comes from transient kinks or bubbles or stems from anharmonic elasticity of DNA requires improved computational methods and further studies." Working in this direction, and being aware of the fact that DNA is in fact an extensible molecule [83], our high cyclization predictions for small lengths in the presence of extension and shear aim to add a piece to the puzzle. Note that this is achieved even under simple linearly elastic assumptions. We believe this mechanism to be relevant and general enough to be shared by several different problems in biology.
Furthermore, birod models [84,85] with sequencedependent parameters are more accurate in capturing DNA conformations, but the theory devised here is comprehensive and can be applied analogously to this level of complexity, allowing the computation of different ring-closure probabilities without involving expensive MC simulations. In the future, in the wide context of end-to-end probabilities, the effect of external loadings will also be investigated. dinates on SE (3). As done originally by Feynman [21], a path integral can be defined via a "time-slicing" procedure, or "parameter-slicing" in our case, which is to replace the infinite-dimensional integral Dq with the limit for n → ∞ of n iterated finite-dimensional integrals n j=1 dq j . These have to be performed on the space of framed curves, whose measure can be chosen to be the product of the Lebesgue measure on the three-dimensional Euclidean space E (3) and of the Haar biinvariant measure on SO (3), which may be uniquely defined up to a constant factor [86,87].
In order to avoid difficulty that can arise from the nonsimple connectivity of SO (3), it is often convenient to consider instead its universal (double) covering SU (2). Any matrix in SU (2) can be parametrized by a quadruple of real numbers γ = (γ 1 , γ 2 , γ 3 , γ 4 ) living on the unit sphere S 3 in R 4 , i.e., γ · γ = 1. The latter quadruple is know as a unit quaternion or a set of Euler parameters [68]. Recalling that by Euler's theorem each element of SO(3) is equivalent to a rotation of an angle ϕ about a unit vector w, the Euler parameters are expressed as a function of ϕ and w as γ 4 = cos (ϕ/2), γ i = w i sin (ϕ/2), i = 1, 2, 3. Hence γ and −γ encode the same rotation matrix, and the correspondence from SU (2) to SO(3) is 2 to 1.
Referring to [41], for parametrizing the group of proper rotations we restrict ourself to one hemisphere of the unit sphere S 3 in R 4 , and we introduce the matrices B 1 , B 2 , and B 3 in R 4×4 : satisfying the algebra B j B k = −δ jk 1 − i jk B i , where i jk is the total antisymmetric or Levi-Civita tensor and summation over equal indices is intended. Furthermore, given a unit quaternionγ, {B 1γ , B 2γ , B 3γ ,γ} is an orthonormal basis of R 4 and each quadruple of Euler parameters γ (hence each rotation) can be expressed in coordinates with respect to the latter basis. In particular, for one hemisphere of S 3 , we consider the new variable b = (b 1 , b 2 , b 3 ) ∈ B 3 1 living in the open ball of R 3 such that γ (b) = 3 i=1 b i B iγ + 1 − b 2γ . Therefore, γ (b) defines a 1-to-1 parametrization of SO(3), adapted to the rotation expressed by the unit quaternionγ, meaning that γ (b = 0) =γ. To be precise, we should remark that the image of such a parametrization does not include the elements lying on a maximal circle (which depends onγ) of the unit sphere in R 4 , since SO (3) is not simply connected and rotations about a generic axis of a fixed angle are inevitably neglected.
For Euler parameters, the infinitesimal measure is given by dq j = δ(1 − γ j 2 )dγ j dr j , so that the Haar volume measure on SO(3) becomes a surface measure on S 3 [86]. Thus, the parametrization φ = γ (b) : B 3 1 ⊆ R 3 → M ⊆ R 4 , with M an hemisphere of S 3 , naturally induces a metric tensor g on the tangent space at each point of M. Denoting the coordinate vectors as φ i = ∂φ ∂b i , i = 1, 2, 3, the components of the metric tensor are given by g i,k = φ i · φ k , i, k = 1, 2, 3, and we get g(b) = 1 + b ⊗ b/(1 − b 2 ), dq j = det [g(b j )] db j dr j with the metric correction being equal to 1/ 1 − b j 2 .
Last, in order to deal with variables defined in the whole of R 3 , we introduce the Gibbs vector c = b/( 1 − b 2 ). As a consequence, we have derived aγ-adapted parametrization of SE (3) denoted by q(s) = (c(s), t(s)) ∈ R 6 as reported in Eq. (14). In particular, exploiting the Feynman discrete interpretation of the path integral measure [21], we obtain Eq. (15).