General theory for plane extensible elastica with arbitrary undeformed shape

A general expression for the strain energy of a homogeneous, isotropic, plane extensible elastica with an arbitrary undeformed configuration is derived. This energy constitutes the correct expression for one-dimensional models of polymers or vesicles, whose natural configuration is characterized by locally changing curvature. We derive the macroscopic stress-strain relations, providing an universal criterion for the neutral curve location. In this respect, we demonstrate that the natural curve existence constitutes the fundamental requirement for the conformational dynamics of any inextensbile biological filament.


Introduction
The old theories of flexure were based on the assumption that the elastica strain consists of extension and contraction of longitudinal filaments, said fibers [1,2,3]. In this context, the reduction of the problem to the one-dimensional plane strain [4] was naturally taken as the compelling point of view. By instance, Lagrange writes about the elastica as a fil flexible et en même temps extensible et contractible [5]. Yet, in Kirchoff theory the three-dimensional deformation of a slender elastic rod is reduced to the bending deformation of a one-dimensional curve [6]. Moreover, the theory of the bending and twisting of thin rods and wires was for a long time developed independently of the general three-dimensional equations of elasticity, by one-dimensional methods akin to those employed by Bernoulli and Euler [7]. However, the problem of which, among the fibers composing the elastica, was eligible to be taken as the representative curve, remained unsettled till the profound works of Parent and Coulomb. Already in 1695, Jacob Bernoulli had argued that the bending moment had to be taken, at each cross-section, with respect to the point where it intersects the line of fulcra, i.e. the line which does not suffer any extension nor contraction all along the deformation: the neutral fiber. If the tensions vary linearly over the rod rigid cross-section, the neutral line coincides with its central line. Although this fact was already known by Beeckman, Hooke, Huygens, Varignon and Mariotte [3,8], it was only in 1713 that these observations took the mathematical form of the Parent criterion [9].
The issue of the existence and ensuing location of the neutral fiber, avoiding any specific elastic hypothesis about the law of variation of the tensions over the cross-section, remains today unjustly unattended by modern scientists. Indeed, what was considered as the elastic problem for a hundred of years has been regarded with condescension by later historians. Mostly because it is wrongly assumed that later works on three-dimensional elasticity somehow fixed it. However, in modern theories of bending, such as Saint Venant's, the existence of the neutral line is implicitly postulated ad hoc, not demonstrated [10,3]. Our purpose, in this paper, is to tackle a problem of a more limited reach: to identify a criterion for the neutral fiber existence and placement, if the linearity of the Hooke law is assumed, but the starting configuration of the elastica is not necessarily straight, as in the slender beam classical theory, but any.
The problem of spontaneously curved bands was already known to Jacob Bernoulli as testified by a posthumous fragment published in 1744 [10]. In the same year also the Euler's treatise on elastic curves appeared, in which the law of an elastica endowed with a non-zero curvature was asserted [3]. A systematic theory of rods with curved undeformed configuration is substantially due to Clebsch [11], although outlined by Kirchoff [6], to whom probably the notion of "unstressed state" must be acknowledged. Only in Timoshenko's theory of curved beam [12], however, it is shown that the neutral fiber does not correspond to the line of centroids. Surprisingly, in the 20th century the problem of the elastica whose undeformed shapes are circular arcs or rings, in particular, has been taken up by several authors [13,14,15,16,17,18,19,20,21,22,23,24], seemingly neglecting Timoshenko's work. As a matter of fact, when the one-dimensional representation of the naturally skew elastica is provided, usually this corresponds to the line of centroids [14,15], assumed, sometimes tacitly, to be the neutral fiber, which can be extensible [14,15,13,16,25,26,27], or inextensible [17,18,19,20,21,22,23,24].
In this paper we consider the deformation of an elastica with a spontaneous curvature which is not constant, thus generalizing the theory of curved beams. Our work is motivated by the fact that natural biological filaments cannot be considered straight in their native state, nor endowed with a uniform curvature. To the contrary, biological materials develop into a variety of complex shapes, which can sense and respond to local curvature. This field of mechanobiology [28] poses fundamental questions about the function of geometry and, particularly, on the role that the local curvature plays in biological systems [29]. In particular, the curvature appears as a determining factor in important biological functional tasks executed by many bio-polymers inside the cell [30,31,32,33,34], as well as in the design of next-generation genomics technologies [35].
On top of that, none of the three major classes of cytoskeletal filaments, microtubules, actin filaments, and intermediate filaments, can be considered homogeneous nor isotropic, as each is formed from chains of discrete protein monomers [36]. At the typical level of description in usual molecular mechanics models, however, they are often treated as traditional engineering elements [37]: they are considered plane elasticae whose mechanical behavior is entirely dictated by the Bernoulli-Euler theory. Our model does not deviate from this conceptual framework; however, our construction can be easily generalized to realistically account for non-homogeneous filaments with varying geometrical and mechanical properties.
Our starting point is the derivation of the strain energy associated to any three-dimensional elastica deformation, as the energy cost of a shape fluctuation is the significant starting concept of any dynamical model in polymer physics [38,39]. The one-dimensional representation of the elastica arises as the natural context for the strain energy formulation (Section 1). The macroscopic constitutive relations are obtained from the strain energy function (Section 2), and in Section 3 we demonstrate that the strain energy minimization furnishes the criterion for locating the neutral curve. We conclude the paper in Section 4, discussing the implication of our findings in terms of polymer dynamics and statitics, providing the evidence that the developed framework yields the correct generalization of flexible and semiflexible polymers models, and of simplified one-dimensional models of fluctuating vesicles. We discuss the relevance that the old-fashioned concept of neutral fiber acquires in modern polymer physics. We show, indeed, that a renovated interest around it, is motivated by the fact that its existence is provided by the minimum energy principle at the base of any polymer conformational change. We discuss the correct and convenient way of implementing our model in computer simulations and, finally, the relevance of our model in describing real biological non-homogeneous filaments.

The Model
We consider the deformation of an homogeneous plane extensible elastica [14]. This is defined as an hypereleastic three-dimensional body of rectangular section, thin in two of its dimensions (height h and depth b), but not in its length, and subject to the only restriction that its deformations take place in a plane, thus excluding twisting. According to the Cosserat's theories [40,41,42,43,44], any elastic deformation stems from an undeformed body configuration or undeformed state, defined as the shape corresponding to the unstressed condition [7]. The work performed in passing from the undeformed to the deformed state is the strain or stored energy, which is the pivotal concept in our model.
We propose to compute the strain energy using a finite difference scheme, as displayed in Fig. 1. The undeformed body is schematically subdivided in N elementary blocks of arbitrary volumes ∆V (i) (i = [1, N ]) (Fig. 1A). To each elementary volume is assigned a local reference system, where the ξ and η α axes define respectively the block's longitudinal and transverse directions, and the origin is arbitrarily placed at an height αh (0 ≤ α ≤ 1) from the block's bottom surface, see Fig. 1A, and at the left side of the block as shown in figure. We now introduce a reference segment of length ∆L (i) α , equivalent to the longitudinal dimension of the plane η α = 0 (Fig. 1A). Assuming the cross-sections to remain rigid, the body deformation necessarily involves the deformation of each block, passing from the volume ∆V (i) to ∆v (i) . At the same time, the reference segment ∆L (i) α transforms into ∆l (i) α . As in the old theories of elastic flexure, we assume that the stress is determined by the strain [45], and the strain consists of the displacements 4 of independent longitudinal filaments u (i) (η α ) = ∆l (i) (η α ) − ∆L (i) (η α ). Embracing the microscopical validity of the Hooke's law, we sum,à la Leibniz, over the contributions of the fibers upon the entire cross section. Thus, the elementary block's strain energy can be written as where Y is a parameter which has the dimensions of a stress and depends on the material properties. Therefore the strain energy (1) represents the work performed in deforming the elementary block from the unstressed configuration reported in Fig. 1A, to that in Fig. 1A ′ . However, in continuity with the old theories, the strain energy (1) is associated to the deformation of the representative material segment. Nevertheless, ∆E α is a scalar quantity and, as such, must be invariant under change of the local reference frame, i.e. change of material representative segment. As a matter of fact, by applying the transformation η α = η β − (α − β)h, the property ∆E β is always satisfied, as it is demonstrated in the Supplementary Informations (SI). The analytical derivation of the energy (1) is fully developed in SI. We hereby report the final expression: where we have introduced the axial strain and the bending measure [14] µ The quantities R (i) α and r (i) α are the radii connecting the intersection point between the extensions of the block limiting sections, with the origin of the block's local reference system: they have a sign corresponding to their orientation, as shown in the SI (see Fig. 1A,A ′ for the details). In the following we will show how, in the continuum limit, they correspond to the local radii of curvature of the undeformed and deformed configurations, respectively. For simplicity, from now on we will refer to such quantities as curvature radii also in the discrete case. Likewise, we will define the block spontaneous curvature as K Most importantly, according to the Timoshenko's curved beam theory [12], F (i) α represents the reduced area, I (i) α the reduced moment of inertia and S (i) α is the reduced axial-bending coupling moment. They have a simple integral expression furnished in SI, corresponding to the formula obtained by Kämmel in the theory of a ring allowing axial compressibility and subjected to hydrostatic pressure [13], and subsequently used in the analysis of compressible rings [15]. Our derivation, however, sheds light on a fundamental aspect that seemingly remained unnoticed in the past analysis: that is, F If the center of the intersection point is placed below the block's bottom line (α = 0), hence K 6 with S α > 0 for α < α U , and S α < 0 for α > α U . Moreover, it is worth to notice that it results lim This means that the fiber allowing the axial-bending uncoupling coincides with the line of centroids in the case of flat undeformed configuration, recovering the classical beam theories result, as shown in the SI. Formally, the full body strain energy is the sum of the elementary blocks strain energy, i.e.
α . In our finite difference scheme, the ordered sequence of consecutive ∆L  α , that assumes an alternative but equivalent expression to that in (4) [14,13,15,46]: where Φ (i) (φ (i) ) is the (i)-th cross-sectional bending angle with respect the x-axis in the undeformed (deformed) configuration (see Fig. 1B and 1B ′ ). Hence In our finite difference scheme, the blocks are assumed to be independent. Thus, nothing prevents us from assigning to each block a different local reference frame, i.e., a different α. As a matter of fact, thanks to the property of invariance of ∆E α under a change of local reference frame, the total energy E α is also invariant. We will revisit this important consideration in Sec. 3 when we discuss the existence of the neutral fiber.
The continuum limit is taken by increasing N up to the point that the poligonal chain L α (l α ) tends to a finite curve, namely the reference curve, or representative fiber L α (ℓ α ) (Fig. 1C, 1C ′ ). According to the older theories of flexure, a material fiber can be envisaged as a translation of one of the outermost curves (α = 0 or α = 1) along the orthogonal cross sections of the body. The Cartesian components of the material line, or fiber, in the laboratory frame are expressed as L α (s) (X α (s), Y α (s)) and l α (s) (x α (s), y α (s)), respectively, where s is the same internal parameter for both the undeformed reference fiber L α : [s m , s M ] → R 2 and deformed one ℓ α : [s m , s M ] → R 2 .
By introducing the tangents T α = dLα ds and t α = dlα ds to L α and ℓ α , the continuum limit of the elastica strain energy is where (·) ′ denotes the derivative with respect to s. The expression (11) costitutes one of the main results of our analysis, extending the Timoshenko's linear theory of curved beams to the case of an elastica with generic undeformed condition, identified by varying local spontaneous curvature K(s) = |Φ ′ (s)|. In SI the differential form of the quantities in (5), (6) and (7) are furnished. In particular, it is worth of noticing how the condition for the axial-bending uncoupling (8) must be valid locally in the continuum limit, i.e. S α (s) = 0. This can be achieved only if α U in (8) changes with s according to The classical case of a thin slender rod, beam or bar [3], is also derived in SI for the sake of completeness. It is shown how this constitutes the zero curvature limit of the general model represented by the energy expression (11).

Macroscopic constitutive relations
We now derive the macroscopic stress-strain relations for the elastica, given a generic undeformed configuration. Our starting point is the strain energy density, defined as W α = ∆Eα ∆Lα . Therefore, when the undeformed state has a constant curvature throughout the whole elastica, from Eq. (2) it turns out that for a generic choice of the representative fiber 8 The macroscopic constitutive equations can be drawn from (13), according to the derivation furnished in [14], which employs the definition of the elastica as a three-dimensional body, or to an alternative variational approach in which the elastica is treated as a one-dimensional medium [25]: N α and M α are the elastica axial force and bending moment respectively, where the one-dimensional representation of the elastica coincides with a generic fiber α. We recall that the bending measure enjoys two analogue definitions, Eq.s (4) and (9). From Eq.s (14), it is immediately verified that linear decoupled relations only hold for the choice of the material fiber allowing for the axial-bending uncoupling of the strain energy (2), namely α = α U in Eq. (8). The same condition is valid in the classical case of flat initial configuration, as it is shown in SI.
It is instructive to see how the constitutive relations transform under change of reference line, i.e. when shifting from the fiber α to β. To this end, we firstly need to express how the strain and the bending measure are modified if the material line changes: and Therefore, plugging these relations into (14), it turns out that The former transformations appear to be universal, in the sense that they are valid in the case of slender bars as well, as shown in SI. Importantly, while the axial force is invariant under change of representative fiber, the bending moment gains an axial force contribution which accounts for the different curvatures of the fibers α and β. Moerover, when N α = 0, the bending moment is invariant under changes of representative fiber: M α = M β . Physically this corresponds to the fact that, for a pure bending transformation, the bending moment cannot depend on the fiber chosen as representative.
For any other transformation for which N α ̸ = 0, the bending moment will depend on the fiber on which it is computed, and the axial force contribution has to be summed up to the total moment of the forces.
When the curvature is locally changing within the elastica, the continuum version of the strain energy density (13) is from which the constitutive equations in a local form can be defined as Once again, the value of α needs to be varied throughout the elastica according to the law (12), for the stress-strain relationships (19) to be locally decoupled.

Neutral fiber location
The neutral fiber is the locus of points which does not undergo any longitudinal extension nor contraction during the elastica deformation [7]. If there is one, this does not necessarily coincide with the line of centroids [7] and its existence depends entirely on the transformation put in use. After the Jacob Bernoulli's first erroneous attempt of establishing a general principle for the neutral fiber placement, based on the balance of bending moments between elongated and compressed fibers, the correct condition was enunciated by Parent in 1713 [9] and rediscovered by Coulomb 60 years later [47]: the neutral fiber is the locus of points separating the region of extension from that of compression, such that the longitudinal forces are balanced.
This criterion, however, is easy to implement in case of Hookean dependence of the tension over the cross-section, but is hard to be assessed for a generic (non-linear) stress-strain relation, as precisely sought by Jacob Bernoulli.
Let us see how Parent's principle applies to the case under investigation. We start by considering an elastica with uniform spontaneous curvature, a situation encompassing the cases of a slender bar, a circular arc or an hyperbola. The requirement of the longitudinal forces balance corresponds to N α = 0. However this constitutes just the necessary condition for the existence of the neutral line. As a matter of fact it does not furnish a precise answer about its placement because, in view of the first of Eq. (17), it does not depend on the specific choice of the representative fiber. However, among all the representative fibers, for rigid cross sections, only one satisfies the zero strain condition ε α = 0, i.e. the value of α U in (8). Therefore, as a general principle, we can conclude that the neutral fiber, when it exists, corresponds to the material representative line which ensures the axial-bending uncoupling of the strain energy. This is an universal property, valid for both flat as well as for non-flat initial configurations. Moreover, it appears that for any transformation guaranteeing the existence of the neutral fiber, the bending moment is independent of the representative fiber choice, as implied by Eq.s (17).
We now investigate how the neutral fiber is identified when the curvature is not constant throughout the elastica in the undeformed configuration. This is achieved by using the macroscopic constitutive relations in the local formulation (19). The condition implied by the Parent principle requires N α (s) = 0 ∀s, accompanied by the zero strain condition |t α (s)| = |T α (s)|. While in the case of elastica endowed with constant spontaneous curvature the neutral fiber is identified with a material line corresponding to a constant α, the situation is totally different for a generic undeformed configuration. As a matter of fact, in this case one cannot properly talk about neutral fiber, in the sense that the curve guaranteeing N α (s) = 0 together with the zero local strain conditions, does not overlap with one of the elastica material fibers. Rather, the concept of neutral fiber is replaced by that of neutral curve, defined as the locus of points identified by the value of α U (s) in Eq. (12). A typical situation is that depicted in Fig. 2C, where the initial curvature is a function of the internal parameter s, K α (s), and the neutral curve is shown as a black dotted line, jiggling around the line of centroids (dashed red line), and approaching the side of the elastica with maximum value of the natural curvature. Notably, while any representative fiber is always orthogonal to the cross-section, the neutral curve is not in general orthogonal.
Therefore, if the elastica transformation provides the existence of a neutral curve, it is possible to adopt a neutral arc-length parametrization such that |T α U (s)| = 1. Adopting this parametrization, the energy assumes the following simple form (see SI) At the same time, the bending moment takes the form Let us take a moment to clarify what is implicit in the expression (20).
The neutral curve requires that the value of α changes with s according to Eq. (12). As a consequence, the neutral arc-length parameterization entails that the strain energy (20) interpolates the expressions (11), in the same fashion as α U (s) interpolates among the different α. However, for a deformation that preserves the neutral curve, the value of E α U ≡ E α , thanks to the property of invariance of ∆E α under a change of reference frame, which holds locally also in the continuum limit.

Discussion
Our one-dimensional treatment of a three-dimensional elastica is in the wake of the old theories of flexure and bending. Nevertheless, the choice of a fiber α as the representative one-dimensional medium is purely arbitrary. As a matter of fact, the linchpin of our analysis has been to disentangle the representative material fiber from the notion of neutral fiber. Indeed, contrary to the modern three-dimensional theories of bending, we do not postulate the existence of a neutral fiber, neither we identify it with the line passing through the centroids. Our main achievement has been to furnish the exact criteria for its existence and its location. These criteria ultimately stem from the derivation of the strain energy, and are encapsulated in the following universal condition: Locating the position of the neutral fiber does not have only an effect on the stiffness of a given cross-sectional form, the Eq. (22) implies that it becomes a crucial issue for any physical transformation aiming at minimizing the amount of work required, given a shape transformation. Let us explain this point with the help of Fig. 2: imagine to have a slender bar, such that in panel A, transformed into the circle in panel A ′ . Among all the possible deformations that bend the elastic bar into a ring, only one allows the presence of the neutral fiber, and is that leaving unvaried the line of the centroids. In this case, Parent's principle is satisfied and the longitudinal forces are balanced, i.e. N α = 0 ∀α. We could imagine to give the bar a final circular shape preserving the length of a different fiber, say the bottom material line (α = 0), together with Jacob Bernoulli (in his first attempt, dating 1694) or Euler [3]. In this case, however, N α ̸ = 0. Therefore the bottom line cannot be identified with the neutral fiber although ε α=0 = 0. The important point is that the amount of work spent by deforming the beam into a circle, is minimum if the line of centroids maintains its length unvaried. This is plainly demonstrated in SI. Now, it is also instructive to consider the opposite case, where the circular shape is the unstressed condition (Fig. 2B), and the flat configuration is the deformed one (Fig. 2B ′ ). The neutral fiber in this case does not correspond to the line of centroids (red dashed line) but it is a circumference whose radius is smaller, according to the value of (8) (black dotted line). Thus the minimum amount of work required to "stretch" the circle into a bar is obtained if the corresponding transformation leaves the neutral fiber unvaried. As a corollary, the bar shown in panel B ′ of Fig. 2 is smaller than that in panel A, because the length of the extended beam in panel So far we addressed the cases of elastica endowed with spontaneous constant curvature. We now turn to the elastica with intrinsic local curvature as in Fig. 2C. The variational version of the minimum principle in (22) is Imagine an elastic filament undergoing a shape transformation. One is generally inclined to think that such a system would prefer to deform by minimizing its energy, skewing and swelling in the fashion that is the easiest of all. The expression (23) specifies that such a minimum-cost deformation is only obtainable by preserving the length of the neutral curve. However, due to the immaterial nature of the neutral curve, its profile may differ considerably from the actual three-dimensional configuration. This is graphically represented in Fig. 2C ′ where, for the sake of simplicity, the elastica deformed configuration is set to be straight and the neutral curve is quite something else.

The importance of the neutral curve for ideal chains
The above considerations have a tremendous impact into both non-equilibrium properties and equilibrium statistics of polymers or other macromolecules. Take for example the conformational dynamics of ideal chains [38,39]. Ideal chains provide simplified one-dimensional models for real flexible and semiflexible biomolecules, being substantially the trace of their backbone, with negligible self-interactions or excluded volume effects. Among semiflexible models, certainly the most successful is the Worm-Like Chain (WLC) introduced by Kratky and Porod [48,49] for inextensible stiff-rod polymers, for which the energy associated with conformational fluctuations may be captured by using merely linear elasticity. The effective WLC energy associated with the bending is [50,51,49] where k B is the Boltzman constant, T is the temperature and A is the persistent length, characterizing the bending stiffness of the polymer. As the expression (24) is derived from thin-rod linear elasticity, it entirely fits into the framework developed in this paper. Indeed the stretching energy of a straight ribbon with inextensible neutral fiber, i.e. the line of centroids, can be expressed as if the same line of centroids (α = 1/2) is assumed as the representative fiber (see SI). The analogy between the expression (24) and (25) is apparent if one recalls that i) both expressions are derived assuming the arclength parametrization for which t (1/2) (s) = 1, ii) the local curvature k(s) ≡ |φ ′ (s)|, and iii) the Young's modulus Y exhibits a temperature dependence [52] (in this respect notice that also A may depend on T [53,54]). At the same time, from the comparison of Eq.s (24) and (25) it emerges that the natural configuration of a polymer is tacitly assumed to be a rigid rod [55,56], not only at T = 0. This contrasts with real polymeric chains, whose microstructural shape is naturally coiled. Taking into account the polymers' natural curvature requires the formal extension of the ideal energy (24) along the line marked by our theoretical analysis. In particular, in view of Eq. (23), the inextensibility requirement attains the higher value of a principle of minimum energy. To be inextensible, however, is none of the polymer material lines, but the neutral curve. This practically translates into fluctuations of the three dimensional polymer conformation and of its contour length, although limited, in contrast with the classical WLC picture. Importantly, it elucidates the contour molecule fluctuations observed in short to moderate length molecules, where L/A ∼ 6 − 20 [57], and it may help to explain the extreme bendability of short DNA [58]. As a matter of fact, the difference between the length of a generic fiber and that of the natural curve becomes more and more apparent for short filaments, or for increasing values of the h/L ratio. Hence, assuming the neutral arc-length parametrization, the correct expression for the bending energy is the Eq. (20), while any other material fiber representation based on constant α requires the general expression (11). This contrasts with flexible polymers for which the concept of neutral curve loses importance and, even assuming the neutral arc-length parametrization (12), the energy is given by

Numerical model for vesicles and polymers
Models for polymers, flexible films, membranes or vesicles often require a numerical implementation, in order to study shapes, fluctuations and dynamics, complementing the statistical mechanics analytical approach [59]. These models in general consist of N impenetrable circular beads connected by links, arranged in a configuration which requires a well defined energy cost. A typical example is furnished by planar thermally fluctuating rings used as simplified models of fluctuating vesicles such as red blood cells. The Liebler-Singh-Fisher (LSF) model [60], by instance, considers planar closed chains including the pressure and a curvature energy term such as where k is the constant bending modulus. Many aspect and variants of the LSF have been discussed in the physics literature using a wealth of analytical and computational techniques [61,62,63,64,65,66,67,68,69], allowing the stretching of the ring [70,71], as well as a spontaneous curvature and locally varying bending modulus [20].
In our theory the strain energy function (11) has been derived by a finite difference scheme, as the continuum limit of the sum of the blocks strain energy (1). Therefore, the discrete form of the Eq. (11) constitutes, without any further approximation, the appropriate expression to be used in numerical simulations (see Fig. 1B,B ′ and SI). However, the fact that the uncoupling value of α U in (8) depends crucially on the local value of the curvature of each block, makes it impossible to adopt a global neutral arc-length parametrization for discrete chains. This means that a discrete decoupled expression for the strain energy function is out of the question, because the discrete chain would result in a discontinuous polygonal, when N is finite (green dashed lines in Fig. 1B,B ′ ). From an energetic point of view, there would be nothing wrong, as the energy would remain the same due to the strain energy's invariance under a change of reference frame. However, from a graphical perspective, it would look awful! On the other side, adopting the constant α description yields the correct and more convenient expression for the numerical implementation of our model, valid for any value of N : We emphasize that this expression holds for any elastic chain, be extensible or unextensible. The unextensibility condition should be enforced by keeping the lengths of the segments at the height α (i) U constants, in any conformational change, and modifying the representative α fiber accordingly. However, the detailed procedure of the numerical implementation of our model will be the subject of an upcoming publication.

Real filaments
Real biological filaments are characterized by local properties that can vary smoothly or abruptly along their contour. Our analysis focused on what we believe to be the most prominent property: the spontaneous curvature. However, geometric and mechanical properties may also characterize different portions of cellular filaments. In reality, filaments can only be considered homogeneous on average, with different heights (h) and depths (b) identifying each discrete protein-based segment. Simultaneously, the intrinsic mechanical properties of each monomer can differ, as in the case of the elastic modulus (Y ), which we assumed to be constant in our analysis.
We believe that the entire framework assembled in this paper can be adapted to handle and accurately describe realistic biological filaments. The reasons lie in two main ingredients of our theory: the discrete scheme and the strain energy invariance under a change of reference frame. The finite difference scheme allows for a more precise description of biological shapes than any continuum theory, such as Cosserats' theory [40,43,44,41,42] because any biological filament, whether it is a microtubule, actin, or intermediate filament, is ultimately composed of discrete units. Since the blocks are independent in our model, the geometrical and mechanical characteristics of biological filaments can be easily implemented locally, while preserving the general form of the strain energy as in (10).
Furthermore, while the concept of material fibers becomes meaningless for non-homogeneous filaments, or at least non-parallel ones, the parametrization in terms of a reference segment α ∈ [0, 1], inherent to any block, will always be possible, giving the concept of fiber a geometrical rather than material meaning. Most importantly, the existence of a neutral curve, allowing the axial-bending decoupled form of the strain energy, must remain valid in any context, highlighting once more the general value of our findings. The elementary block, with its integral system of reference ξ-0-η α , has the origin placed at an height αh from the bottom surface (α ∈ [0, 1]). A reference fiber is characterized by its length (η α = 0) is ∆L α (dotted black line) and by R α , with the corresponding spontaneous curvature K α = 1 |Rα| . R α is the distance of the origin of the local reference from the intersection between the extended sides of the block, with a positive sign if its orientation is concordant with that of η α . It becomes the curvature radius in the continuum limit. A ′ : The block undergoes a deformation to a new shape (dark blue), where the original configuration is in light blue. Shape change entails R α → r α . Any fiber of the block endures a deformation represented by the fiber displacement u(η α ) = ∆l(η α ) − ∆L(η α ) (green arrows). B: The undeformed elastica is displayed as a sequence of elementary identical blocks, each with its own shape. The sequence of the reference segments is shown as a dotted black polygonal L α , while bending angle Φ (i) of each block is in red. The green dashed lines correspond to the α (i) U fibers. Conventionally, the construction assumes that the left section of each block remains straight-angled, so that the bending angles are defined at the left side of each block. B ′ : The deformation of the elastica involves the deformation of the reference polygonal chain, L α → l α (dotted black line). The bending angles vary from Φ (i) to φ (i) . C: In the continuum limit, the material reference curve is represented by a solid black curve. The tangent to this curve, T α (s), forms with the x axis the spontaneous bending angle Φ(s). C ′ : The deformed material curve (solid black line) stems from its undeformed configuration (dotted black line). The black arrow represents the new tangent t α (s), while its corresponding bending angle φ(s) is in red. The beam in panel A is bent into a ring. If the line of centroids maintains constant its length, therefore this coincides with the neutral fiber and the axial forces are balanced (N α = 0, ∀α). Among all possible deformations which transform a beam into a circle, the one allowing the existence of the neutral fiber costs the minimum amount of work. B: The undeformed condition is a circle. The line of centroids (α = 1/2) is depicted as a dashed red circle, while the material fiber granting the axial-bending uncoupling, defined by α U in Eq. (8), is shown as a dotted black line. B ′ : The circle in panel B is deformed into a straight beam. If the fiber identified by α U keeps its length constant during the circle extension, therefore it coincides with the neutral fiber and the Parent's principle is satisfied. Any other stretching of the circle into a bar will cost more energy than that leaving unchanged the neutral fiber. Notice that, although the circles in panels A ′ and B have the same dimensions, the absence of axial forces implies that the beam in A is longer than that in B ′ (see SI). C: The generic undeformed configuration has a line of centroids shown as a dashed red line. The curve allowing the axial-bending uncoupling in the Eq.(11) is show as black dotted line. This curve does not overlap with any material fiber according to (12) (1), where the quantity Λ(ηα) are represented by orizontal thin red arrows. The spontaneous radius of curvature relative to the representative fiber Rα is shown by a vertical thick red arrow. A ′ : In the deformed condition, the size of the representative fiber varies from ∆Lα to ∆lα and its radius of curvature becomes rα. In the reference frame integral with the block, any fiber undergoes a displacement u(ηα) (green arrows). At the same time, the size of any deformed material fiber can be expressed as ∆l(ηα) = ∆lα + λα(ηα) (red arrows).
We consider the elementary block depicted in Fig. 1. The undeformed condition is represented in panel A, whereas its deformed state is shown in panel A ′ . The length of the representative fiber in the natural state is ∆L α , and that of a generic fiber at an height η α is The corresponding deformed quantities are ∆l α and ∆l(η α ) are defined by Hence, the deformation that a generic material segment experiences passing from the state A to the state A ′ is expressed as The geometry of the undeformed shape satisfies the following equality [1,2] while for the deformed configuration we have Thanks to Eqs. (4) and (5), the displacement (3) attains the final form By substitution of Eqs. (1) and (6) into we obtain the expression: where (10) and These three factors have different functional forms according to whether R α is positive or negative (see Fig. 2). The quantity R α can be positive or negative, as required by the consistency of the Eq.(4). |R α | becomes the radius of the osculating circle which locally approximates the reference segment in the continuum limit, and the sign of R α is assigned in the following way. It is clear that the intersection C between the sidelines containing the block's sections lies on the axis ξ = 0 of the local reference system ξ-0-η α . If C lies below the bottom fiber α = 0, then R α is positive ( Fig.2A). If conversely C is above the upper fiber α = 1, then R α is negative ( Fig.2A ′ ). It follows that C has coordinates (0, −R α ) in the local reference system. The same prescription for the sign applies to r α . If R α > 0 the solutions of (9), (10) and (11) read and Figure 2. Sign of the radius of curvature. The radius of curvature has a sign assigned whether its direction, connecting the intersection C of the sidelines containing the block's sections with the beginning of the reference segment, coincides with that of the ηα axis of the frame integral with the block (panel A), or opposite to it (panel A ′ ).
with R α = R 0 + αh. When we consider the case R α < 0, the three integrals 10-11 can be solved yielding and where The two opposite bending states in Fig.2 are expressible in the compact forms provided in the main text recalling that, if R α > 0, hence necessarily R 0 < R 1 (K 0 > K 1 ), whilst, for R α < 0, therefore |R 1 | < |R 0 | and K 1 > K 0 . We stress the fact that the expressions of F α , S α and I α are fully established once the value of α and max[K 0 , K 1 ] are furnished. In particular the following identity holds if K 0 > K 1 while for K 1 > K 0 shown as green dots. The spontaneous radius of curvature can be positive (as R 1 α ) or negative (as R 2 α ). The connection between the size of representative segment, ∆L (i) α , and its spontaneous radius of curvature R (i) α , is given by the relation (23), where the bending angles Φ (i) define the deflection of the undeformed blocks from the external x axis. A ′ : each block composing the beam undergoes a deformation, such that the deformed polygonal chains lα is constructed by the sequence of the representative segments whose vertices are x According to the integral expressions (9) and (11), it is always F α > 0 and I α > 0 because the integrand functions are strictly positive in the integration interval. To the contrary, the sign of S α can vary according to the value of α and to max[K 0 , K 1 ] (min[R 0 , R 1 ]). For the sake of clarity, S α > 0 for α < α U , and S α < 0 for α > α U . The choice of α U ensuring the stretching-bending uncoupling, guarantees that and The finite difference scheme, outlined so far, requires the evaluation of the discrete strain energy α needed to deform the elastica in Fig. 3 from A to A ′ . Moving to the laboratory frame we find that the relation is always satisfied [3], where ∆Φ (i) = Φ (i) − Φ (i−1) and Φ (i) is the i-th cross-sectional bending angle with respect the x axis. Since we assume the limit of small deflections, we can approximate tan ∆Φ (i) ≃ ∆Φ (i) . The reference undeformed polygonal chain L α is specified by the series of points L Fig.3A).
In the deformed state (Fig.3A ′ ) we have where ∆φ (i) = φ (i) − φ (i−1) , and φ (i) corresponds to the bending angle between the i−th block and the x axis. Again we assume the small deflection limit, i.e. tan ∆φ (i) ≃ ∆φ (i) . The deformed chain l α has the points We introduce the strain measure as ε while the bending strain measure can be obtained by two different definitions: the first is due to Kammel [4] and the second to Antman [5] One can easily see that they are equivalent by the construction in Fig.6. Upon summation of the terms in Eq. (8), using the definition (26), we obtain the energy The expression of F α , S α and I α in terms of ∆Φ, ∆L α and α, is achieved by inserting the relation (23) into the expressions (12)- (14) and (15)- (17): and We recall that it is convenient to take ∆L α = ∆L 0 −αh∆Φ if min[∆L 0 , ∆L 1 ] = ∆L 0 , and ∆L α = ∆L 1 +(1−α)h∆Φ if min[∆L 0 , ∆L 1 ] = ∆L 1 . The differential strain energy E α is derived by firstly introducing two parametric expressions for the undeformed and deformed reference material curves as L α : [s m , s M ] → R 2 and ℓ α : [s m , s M ] → R 2 respectively. Secondly we take an arbitrary partition s m = s 0 < s 1 < s 2 < · · · < s N = s M to which we connect the two polygonal chains L α and l a , in such a way that the vertices satisfy L α (s i ) ≡ L In the continuum limit, N is increased until the lengths of the polygonal chains L α and l α equal those of the curves L α and ℓ α . This condition is mathematically enforced by the limiting relations as ∆s i → 0. Correspondingly, the two tangents to the curves are defined as T α (s) = dLα ds and t α (s) = dlα ds . Yet, the limit ∆s i → 0 entails The differential strain measures follow from the limit of Eq.s (31) and (32): The differential formula for the three factors F α (s), S α (s) and I α (s) are obtainable from the Eq.s (28), (29) and (30): and which is also expressible as The two formulae (41) and (42) are equivalent and can be obtained recalling that |T 1 (s)| = |T 0 (s)| − hΦ ′ (s). Therefore the uncoupling condition is a local property of the elastica and, for either choice of α U (s), the three factors (38)-(40) reduce to S αU (s) = 0 (44) and The neutral arc-length parametrization requires that |T αU (s)| = 1.
Therefore, adopting this parametrization and the neutral curve as representative of the whole elastica we have that for a generic transformation the energy is expressible as In red is represented the reference frame ξ-0-ηα integral with the block, when the reference material segment is placed at an height αh from the block's bottom surface. The length of the reference segment is ∆Lα. When the reference material segment is placed at a different height βh, the corresponding frame ξ-0-ηβ is depicted in green and its size is ∆Lβ. A ′ : The block undergoes a deformation from its natural shape (light blue): the size of the reference fiber changes to ∆lα or ∆l β . The energy cost associated to this deformation is the same whether the reference segment is placed at αh or βh.

II. ENERGY INVARIANCE UNDER CHANGE OF REFERENCE FRAME
Consider the Fig. 4. The undeformed elementary block is represented on the left A, and its shape upon deformation is displayed on the right A ′ . Let us assign to the block an integral planar reference system, where the ξ and η α axes define respectively the block's longitudinal and transverse directions. The origin of such reference system is placed at an height αh (0 ≤ α ≤ 1) from the block's bottom surface, and on the left lateral block boundary. The quantity ∆L(η α ) corresponds to the length of a generic material line placed at the height η α , with −αh ≤ η α ≤ (1 − α)h. By definition, the value of ∆L(η α = 0) is the representative material line length ∆L α . Now, a translation of the reference system along the axis ξ = 0 is equivalent to a linear change of variables: where η β is the new axis pointing along the block transverse direction. However, the material line lengths ∆L(η α ) have to be invariant under the transformation (49): It is also clear that in the reference system ξ-0-η β , the representative material line has length ∆L β = ∆L (η β = 0) = ∆L (η α = (β − α)h) and −βh ≤ η β ≤ (1 − β)h. Let us turn to the deformed configuration A ′ . The deformed longitudinal length ∆l follows the same law (50) under the shift of the reference system, i.e.
Therefore, since the extension is defined as u = ∆l − ∆L, thanks to Eqs. (50) and (51) we have that the following equality holds The strain energy of the block calculated in the ξ-0-η α reference frame is By applying the change of variables (49), the equality of the integral expression (53) ∆E α = ∆E β follows from the relations (50) and (52). Let us consider the deformation of the elementary block presented in Fig.5. The undeformed flat condition is depicted by a dotted black line, and it has the peculiarity that the longitudinal length is equal to ∆L for any choice of the representative segment. The plane integral reference frame is identified by the ξ and η α axes, pointing respectively towards the block's longitudinal and transverse directions. The origin is placed at an height αh (0 ≤ α ≤ 1) from the block's bottom surface, and on the left lateral block boundary. Any fiber placed at an height η α attains a length ∆l(η α ) upon deformation, with ∆l(η α = 0) = ∆l α . According to the geometrical construction in Fig.5, the Eq.(2), is equivalent to ∆l(η α ) = ∆L + u(η α ).

III. DERIVATION OF THE STRAIN ENERGY FROM A FLAT UNDEFORMED CONFIGURATION
From Eqs. (54), (2) and (5) the elongation of any fiber can be expressed as In this condition, the strain energy of the block takes the form where we have inserted the relation (55). Solving the integral and defining the strain as ε α = ∆lα−∆L ∆L , we arrive at the expression It is clear that the value of α which guarantees the axial-bending uncoupling is the line of the centroids α U = 1/2. The strain energy of the whole elastica is given by the sum over the block contributions, formally by We aim at furnishing, however, its analytical expression in the laboratory frame (Fig. 6). The single block's representative segment size ∆l (i) α is a positive quantity, being ∆l are the vertices of the reference polygonal curve l α in the lab reference system (Fig. 6). The polygonal curve is defined as the ordered sequence of the representative segments ∆l (i) α . The quantity r (i) α , on the other side, can be positive or negative. This is required by the consistency of the Eq.(5). As a matter of fact, r (i) α becomes the radius of the osculating circle which locally approximates the reference segment in the continuum limit, and the sign of r (i) α is assigned in the following way. It is clear that the intersection between the sidelines containing the block's sections lies on the axis ξ = 0 of the local reference system integral with any block (see Fig. 2). Hence, in this reference system the coordinates of the circle's center are defined as 0, −r (i) α : this establishes uniquely the sign of r (i) α . Now, moving to the laboratory frame we find that the relation (24) is always satisfied, with ∆φ (i) = φ (i) − φ (i−1) , and φ (i) corresponds to the bending angle between the i−th block and the x axis. In the small deflection limit, i.e. tan ∆φ (i) ≃ ∆φ (i) , the discrete strain energy for the entire slender beam is therefore framed as The bending measure is defined as or, thanks to (24), as Using these definitions, the energy (61) takes the following form Let us introduce two parametric expressions for the undeformed and deformed reference plane curves as L α : [0, L] → R 2 and ℓ α : [0, L] → R 2 respectively. As a consequence, the Cartesian coordinates of the undeformed reference curve in the laboratory frame are and those of the deformed curve ℓ α (s) are l α (s) ≡ (x α (s), y α (s)). Taking an uniform partition of [0, L], i.e. 0 = s 0 < s 1 < s 2 < · · · < s N = L such that s i − s i−1 = ∆s ≡ ∆L for any i, we obtain that the polygonal vertices are l α (s i ) = l (i) α and the local longitudinal strain ε (i) α is given by Analogously, the bending angles at the polygonal vertices are φ(s i ) = φ (i) . The continuum limit is taken by increasing N until the length of the polygonal chain l α approaches from below that of the curve ℓ α , i.e. |∆lα(si)| ∆s → |ℓ ′ α (s)| as ∆s → 0. The curve derivative is defined as ℓ ′ α (s) ≡ t α (s), where we have introduced the tangent of the curve t α (s) = dlα(s) ds . Finally, if the continuum limit entails that ∆φ(si) ∆s → φ ′ (s) and

IV. STRAIN ENERGY LIMITING CASES: REGAINING THE FLAT UNDEFORMED CONDITION
In the present section we show how to recover the straight beam strain energy (57), from the energy (8) calculated from a generic undeformed configuration. To this aim, it will be sufficient to study the behaviour of F α , S α and I α in the limit of h min[|R0|,|R1|] → 0. Let us firstly express the relations (9)-(11) as and Then we consider the condition h min[|R0|,|R1|] ≪ 1 and expand the logarithm to the third order: Hence we get Hence by subtracting the expression for M β from (76) and recalling the axial force invariance we have The case of a general undeformed condition can be determined as follows. First we notice how the strain transforms under change of material line Moreover, the reduced area, the reduced axial-bending coupling moments and the reduced moment of inertia change under material line transformation as Inserting the previous relations into it easily turns out that the equality N α = N β holds also in this case. Moreover, plugging the same transformations into the second of Eqs.(82), we recover the Eq. (77) also in the case of generic initial conditions.
Without loss of generality, let us adopt the line of centroid as the representative material line, namely α = 1/2. We know that this choice has the only advantage of yielding the axial-bending uncoupling in Eq.(87): Nonetheless, if the transformation is such that r 1/2 = L 2π , i.e. the middle fiber maintains its length constant (zero strain condition), the energy has a minimum. In other words, among all the possible deformations that transform a bar into a circle, that one which leaves unvaried the middle fiber (the neutral fiber) costs the minimum amount of work: This minimum principle can be seen as the straightforward application of Parent's principle N 1/2 = 0. Now let us consider the opposite situation, where a naturally curved beam is flattened into a bar as in Fig.7B-B ′ . The undeformed configuration is given by with θ ∈ [0, 2π) being the internal parameter which is now adimensional, rather than having the dimension of an internal length.
so that Φ(θ) = π 2 − θ. On the other side the equation for the deformed bar of length l is and φ(θ) = 0. Hence, the energy cost connected to such a transformation is In analogy to the previous case, we choose the value of α which entails the axial-bending uncoupling, namely, according to Eq. (41), Thanks to the fact that R α = R 0 + αh, from (95)  Again, the minimum of the energy is consistently required by the validity of Parent's principle. So far we have considered the generic situation where circle in panel B of Fig.7 and that in panel A ′ are different. If we set the same dimensions for both of them, we have R 0 = L 2π − h 2 that, inserted into the energy expression (96) yields E αU (L; l) = πbY l 2 4π 4 ln The amount of work needed to stretch the ring out, keeping constant the length of the neutral fiber, is Conversely if we want to stretch the ring keeping constant the line of the centroids, therefore it is sufficient to replace l = L into the expression (98):