Mechanics of tubular helical assemblies: ensemble response to axial compression and extension

Nature and technology often adopt structures that can be described as tubular helical assemblies. However, the role and mechanisms of these structures remain elusive. In this paper, we study the mechanical response under compression and extension of a tubular assembly composed of 8 helical Kirchhoff rods, arranged in pairs with opposite chirality and connected by pin joints, both analytically and numerically. We first focus on compression and find that, whereas a single helical rod would buckle, the rods of the assembly deform coherently as stable helical shapes wound around a common axis. Moreover, we investigate the response of the assembly under different boundary conditions, highlighting the emergence of a central region where rods remain circular helices. Secondly, we study the effects of different hypotheses on the elastic properties of rods, i.e., stress-free rods when straight versus when circular helices, Kirchhoff’s rod model versus Sadowsky’s ribbon model. Summing up, our findings highlight the key role of mutual interactions in generating a stable ensemble response that preserves the helical shape of the individual rods, as well as some interesting features, and they shed some light on the reasons why helical shapes in tubular assemblies are so common and persistent in nature and technology. We study the mechanical response under compression/extension of an assembly composed of 8 helical rods, pin-jointed and arranged in pairs with opposite chirality. In compression we find that, whereas a single rod buckles (a), the rods of the assembly deform as stable helical shapes (b). We investigate the effect of different boundary conditions and elastic properties on the mechanical response, and find that the deformed geometries exhibit a common central region where rods remain circular helices. Our findings highlight the key role of mutual interactions in the ensemble response and shed some light on the reasons why tubular helical assemblies are so common and persistent.


Introduction
Many biological structures can be modeled as tubular assemblies of helical rods, such as the tail sheaths of bacteriophage viruses [1,2], the cellulose filaments in the tendrils of climbing plants [3], the bundles of microtubules and motors in all eukaryotic flagella and cilia [4,5], the envelopes of shapeshifting unicellular organisms such as Lacrymaria Olor [6] and the pellicle of euglenids, a family of unicellular algae [7][8][9][10] (see Fig. 1).
From a technological viewpoint, tubular structures composed of helical fibers exhibit programmable shape-shifting capabilities. This makes them adaptable to changing functional needs, by varying their conformation and properties. Exploiting these features, they have been employed in a broad range of domains, e.g., deployable antennas in aerospace engineering [11], tubular vascular stents in biomedical engineering, sheaths of McKibben artificial muscles and, more generally, helically-arranged fibers in soft robotics and biorobotics [12][13][14][15]. In all these examples, the system consists of a tubular network of helical rod-like structures. Understanding the mechanical behavior of these systems requires knowledge of how single helices deform under external loads, and how they interact with each other to produce the ensemble response.
In this work, we address these questions by focusing on the response under compression and extension of a tubular assembly made of 8 helical rods. The specific geometry of this assembly, inspired by Ref. [11], is further described below. The behavior of the assembly is analyzed through numerical simulations, adopting a Kirchhoff rod model, that is, an inextensible and unshearable Cosserat rod [16]. For this purpose, the finite element software COMSOL Multiphysics ® v5.4 is used in equation mode. Specifically, the principle of virtual work is implemented to model the nonlinear response of interacting helical rods in the large deformation regime, under prescribed end-displacements. Fig. 1 a Neck protrusion with shape changes in the microtubule meshwork of Lacrymaria Olor [6]. b Metaboly of Euglena Gracilis [10], a shape changing mechanism based on sliding helices. Adapted with permission from Ref. [8] Helical structures have been studied extensively and, in particular, their treatment using Kirchhoff rod theory has been discussed in many textbooks and research papers (e.g., see Refs. [16] and [17] and the many references cited therein). Recent contributions to the literature on the mechanics of helical rods include [18][19][20][21][22]. With respect to this recent literature, the main differences of our research consist in the type of structure analyzed, and in the generality of the allowed deformations, i.e., rods are not assumed to remain circular helices a priori and can deform into helical shapes with nonconstant curvature and torsion. Another relevant source is the book by Costello [23], where the behavior of wire ropes made of bundles of helical strands is investigated. The main novelties introduced in our work are the type of assembly considered, the presence of large deformations, and the study of the effect of boundary conditions in structures of finite length, whose deformed shapes may deviate from that of circular helices.
Concretely, in this paper we consider a tubular assembly consisting of 8 rods, arranged in 4 pairs of circular helices with opposite chirality and same axis, connected by pinjoints. Friction between the rods, contact and deformations of the cross sections are not accounted for. Rods are subjected to boundary constraints only at their edges. We first examine the response of the assembly in compression and extension in some specific examples that highlight how the rods in an assembly deform all in the same fashion, i.e., as a coherent ensemble. Then, focusing on the behavior in compression, we analyze how the geometry of the deformed configurations is affected by different constraints. We find that, depending on boundary constraints, the rods can either remain circular helices along the whole length of the assembly, or develop edge layers with perturbed geometry around a bulk region where circular helical shapes persist. We then study how the relation between axial force and axial strain is affected by different hypotheses on the elastic properties of the individual rods. In particular, we consider rods that are stress-free in the initial configuration (as circular helices) and rods that are stress-free when straight. Finally, we analyze the effect of modeling the rods as ribbons, and find that this may produce dramatic changes in the stiffness of the assembly, which become evident in a numerical extension experiment.
Our results show that, when arranged in tubular assemblies, helical rods stabilize each other, and that they respond collectively in a way that may be very different from the response they would exhibit if they were loaded individually. This stable collective response is made possible by the mutual interactions between rods, represented by the reactions at pin-joints. This is particularly evident in compression experiments: an individual helical rod typically buckles away from the helix axis; on the contrary, in the assembly, rods deform in such a way that the structure remains symmetric with respect to its initial axis. The deformed shapes of rods with same (opposite) chirality are then identical modulo rigid rotations (reflections). In other words, the rods of the assembly deform cooperatively into helical structures wound around one common axis, that is, the initial axis of the assembly. This collective behavior promotes the stability and persistence of helical shapes in a tubular assembly and, in turn, suggests why these shapes are so commonly observed in nature and technology.
The remainder of this paper is organized as follows. Section 2 introduces the system under investigation and describes its model in terms of kinematics and equilibrium equations, and it delves into the theory of cylindrical deformations; it then describes the finite element method (FEM) implementation of our model and it outlines the investigated problems where different boundary conditions are applied. In Sect. 3, we discuss the main results of our study: we compare the solutions for several boundary conditions, and we analyze cylindrical deformations in the case of zero intrinsic curvatures and torsion, with an example of a zero-stiffness structure. Moreover, the impact of choosing an alternative constitutive model for the rods, namely Sadowsky's ribbon model, is examined. Section 4 reports our conclusions and outlook for future work.

Analytical formulation and computational model
We consider an assembly of initial height h 0 composed of 8 helical rods, i.e., they can be described as circular helices when unloaded. Our goal is to study its equilibrium configurations under compression or extension, enforced by prescribing its height h, with no other forces applied. The 8 helices are made of the same material, they have the same cross-section and length L, initial radius a 0 , initial pitch p 0 , and same axis, as shown in Fig. 2. The rods are organized in 4 pairs of 2 helices with opposite chirality, that is, a right-handed and a left-handed one. In each pair, the terminal ends (tip and base) of the two helices are linked together, and they initially lie on the circumference given by the initial radius. Terminal ends of different pairs are distinct. At crossings, helices are linked to each other by pin-joints, which allow relative rotations along the cross-sectional axis with the largest moment of inertia. Fig. 2 Schematics of the initial configuration of the assembly. The structure is composed of 8 inextensible helices, 4 right-handed (blue) and 4 left-handed (red), connected by pin-joints (green dots). All helices are made of the same material and have the same length, radius a 0 , and pitch p 0

Geometry and kinematics
Let {e 1 , e 2 , e 3 } be the global Cartesian reference frame, with e 3 along the axis of the assembly. Following the Kirchhoff rod theory, each rod is described by a curve and a local three-dimensional reference frame along it. For all helices, we choose the curves joining the centroids of their crosssections as reference curves. In their initial configurations, a right-handed and a left-handed helix belonging to the same pair can be described by the curves c 0 = a 2 0 + b 2 0 , and θ 1 , θ 2 are offset angles with respect to e 1 . The 8 helices composing the assembly have the following offset angles: We can then define a local orthonormal reference frame for each helix, given by where e R r (s) and e R φ (s) (e L r (s) and e L φ (s)) are the radial and tangential unit vectors for the right-handed helix (left-handed helix) at arc-length s, respectively. Using this notation, Eq. (1) can be simplified as We then define the directors of the Kirchhoff rod, choosing d R,L 3 (s; h) equal to the tangent unit vector to the solution curves γ R,L (s; h) and with the generic position r(s; h) recon-structed as In the initial configuration, directors are described by and we can compute initial curvatures and torsion as

Parametrization in terms of quaternions
In numerical simulations, we parametrize directors via quaternions [24], i.e., q(s) = (q 1 (s), q 2 (s), q 3 (s)), q 4 (s), which are four redundant parameters used to avoid singularities. Starting from the identity the rotation tensor R can be expressed, according to Rodrigues's formula, as where a and α are the axis and angle of rotation at s, respectively, and A is the skew-symmetric tensor corresponding to a. Quaternions are then defined as and are therefore constrained by which causes the dependency of one component from the others. We can then rewrite Eq. (10) as where Q is the skew-symmetric tensor corresponding to q, and we finally express directors as

Balance of forces and moments
In the absence of external actions, the balances of forces and moments read where M is the internal moment and T is the internal force. The internal moment for a pair of rods (one right-handed and one left-handed) can be written as We assume the following constitutive equations: are the intrinsic curvatures and torsion of each right-handed (left-handed) helix, i.e., those exhibited in the rest configuration; B 1 , B 2 represent the bending stiffnesses of each rod along d R,L 1 , d R,L 2 , and B 3 their torsional stiffness along d R,L 3 . According to linear elasticity, we have that where E Y represents the Young modulus of each rod, I 1 , I 2 the principal moments of inertia of its section along d R,L 1 , d R,L 2 , G its shear modulus and J its torsional constant. Since the initial configuration is stress-free, we substitute the intrinsic curvatures and torsion with the ones calculated in Eq. (7), that is,

Theory of cylindrical deformations
Given a helical assembly, cylindrical deformations are families of configurations parametrized by the axial strain ε = h/h 0 −1, where each rod behaves as a circular helix, i.e., a rod with constant non-zero torsion, one constant non-zero curvature, and the other one vanishing. For such deformations, the configurations of a rod can be described by the maps in so that the energy of each helix reads which is independent of chirality, because of Eq. (7). The length of each rod can be expressed as with n the number of turns of the helix, which is constant because of the symmetry of the problem. Then and because of the inextensibility, we obtain Rearranging Eq. (4), we can express b for both right-handed and left-handed helices as and exploiting Eq. (26), we obtain According to the Kirchhoff rod theory and Eqs. (27) and (28), we can write For what concerns the loads, substituting Eqs. (19) and (29) into Eq. (18), and considering Eq. (22), we get The internal moment at s = L is equal to the external one applied on the considered helix, and the sum of their internal moments at s = L corresponds to the external moment Q(h) applied to the pair, that is, Since the rods deform as circular helices, the directors are given by Eq. (6), with h replacing h 0 . Moreover, due to symmetry we have that so that we can rewrite Q(h) into the local reference frame {e r (L), e φ (L), e z } as Under the assumptions of circular helices and of a symmetric structure, only radial rotations are allowed at terminal ends, so that the external work performed by moments vanishes: Hence, the only external load expending work is the external force needed to balance the reaction to height variations. From the energy balance of the assembly the external force F(h) tot can be computed as The energy of the system can be written as the one of a single helix times the number of helices, regardless of chirality. Therefore, the average force per helix F z is computed by differentiating Eq. (23), that is, We differentiate curvature and torsion as We substitute Eqs. (29) and (40) into Eq. (39), which reads Finally, according tô we non-dimensionalize Eqs.
where K 3 = B 3 /B 1 . By setting the intrinsic curvatures equal to the initial ones, according to Eq. (47) we haveF

Weak-form of the balance equations
In order to find the equilibrium configuration of the assembly at a prescribed height, we enforce the balance of forces and moments according to the weak formulation of Eq. (16) given by where δ E is the virtual variation of the internal energy of the system, andq 1 ,q 2 ,q 3 ,q 4 are test functions for the quaternions. For each rod, using constitutive equations we express its elastic energy as and applying the virtual variation operator δ we obtain ,Ω (i) are the curvatures and torsion corresponding to the test functions of the quaternions. The additional constraint of Eq. (11) can be imposed using a Lagrange multiplier η (i) . Hence, the total energy of the system reads with virtual variation resulting in Thus, the final equilibrium equation reads

Constraints and boundary conditions
We investigated the response of the assembly as we vary the height h of the assembly by imposing the e 3 -component of the position vector at the terminal ends of the helices. We studied four cases, corresponding to the following sets of boundary conditions applied at the edges s = 0, L of each rod: 1. displacement components along e 1 , e 2 are set to zero, rotations are set to zero; 2. displacement components along e 1 , e 2 are set to zero, only radial rotations (i.e., with axis e R r = e L r ) are allowed and free; 3. displacement components along e 1 , e 2 are free, only radial rotations are allowed and free; 4. displacement components along e 1 , e 2 are set to zero, rotations are free.
Rotations are prescribed by setting Dirichlet boundary conditions on the quadruplet q 1 , q 2 , q 3 , q 4 , specifying appropriate constants. The components of the position vector are imposed by choosing integration constants to reconstruct r(s; h) from the tangent vector, see Eq. (5), or by enforcing weak constraints using Lagrange multipliers, i.e., by adding the term to Eq. (55), where λ i and r i are the Lagrange multiplier and the imposed value for i-th coordinate of r(L), respectively. In particular,r(L) is expressed in terms of (q,q 4 ) by applying the virtual operator δ to Eq. (5). Since the resulting equation must hold for every λ i , the constraint equation r i (L) = r i is enforced for every i. In all cases we imposed no relative translation between helices at the pin-joints connecting them. Moreover, since directors d R 2 , d L 2 must remain aligned at crossings, relative rotations were allowed only about the axis of the pin-joint. These constraints were enforced using Langrange multipliers, by imposing zero relative displacement between helices at pin-joints and the orthogonality of d R 2 with respect to d L 1 , d L 3 , which corresponds to the parallelism of d R 2 , d L 2 .

Numerical aspects
The governing equations of the rod model, along with boundary conditions, were implemented and solved in COMSOL Multiphysics ® v5.4. In particular, the weak form PDE mode was adopted, using a custom implementation of the equations of the model. Regarding the mesh, upon checking that the numerical solution was insensitive to further refinement, each rod of the assembly was divided into 10 elements. Quadratic shape functions were used to discretize the quaternions fields.
In all the boundary value problems a continuation approach was followed, so that the position of one boundary point was varied gradually. Indeed, at each boundary displacement step, the system of non-linear algebraic equations resulting from the finite elements procedure was solved by a quasi-Newton iterative algorithm, taking as initial guess the numerical solution at the previous step. A direct solver (PARDISO) was chosen to solve the linearized system obtained at each continuation step.

Results and discussion
We report and discuss the main findings on the response under compression and extension of the assembly described in Sect. 2, as extracted from FEM simulations. We compare computational results with analytical predictions, and we explore the special case of naturally-straight rods, i.e., with κ s 1 = κ s 2 = Ω s = 0, and the consequences of a different constitutive model, that is, the Sadowsky's ribbon model. To allow meaningful comparisons, the main variables are non-dimensionalized according to Eqs. (43)-(45).
Inter-penetration is not explicitly forbidden in our model, and this is clearly an issue in regimes of extreme compression. Indeed, our ideal assembly is allowed to degenerate into a single circular ring while, in the physical system, the thickness of the rods would put a limit to the minimal height to which the ensemble can be compressed. To thoroughly explore these regimes, our model should be extended to take self-contact into account.

Stabilization effect due to ensemble response
Single helical rods are known to be unstable under compression, with critical loads determined by their slenderness and material properties. It is then remarkable how helices, when they form assemblies, tend to stabilize each other through mutual interactions, which collectively determine the ensemble response of the structure. Indeed, the higher stability of helical assemblies is arguably a key factor for their recurrence in natural and artificial structures. In this regard, Fig. 3 compares the behavior of a single right-handed helix and of the 8-helix assembly at ε = −0.5, with boundary conditions corresponding to Case 1. The single helix buckles away from the vertical axis, that is, its helix axis in the initial configuration; on the contrary, in the assembly, rods deform in such a way that the structure remains symmetric with respect to its initial axis. Therefore, because of axisymmetry, the deformed configuration of rods with same (opposite) chirality can be obtained from one another via a rotation (reflection). Thus, the mid-line of each rod remains confined to a surface of revolution around the initial axis of the assembly. The resulting envelope surface, termed in this way because it represents an envelope for the deformed shape of all rods, provides a useful tool to describe the way in which rods deform cooperatively and collectively, i.e., as an ensemble, into helical structures wound around a common axis, given by the initial axis of the assembly.
This type of collective mode of response, in which the deformed configurations of the individual rods are identical modulo rigid rotations or reflections, arises also under extension, as reported in Figs. 4 and 5. In the next subsection, we examine in detail the impact of different boundary conditions on the response of the assembly in the case of compression, in which we can observe the most dramatic differences between the response of a single helix and the one of the assembly (compare Figs. 3 and 4).

Comparison of equilibrium configurations corresponding to different boundary conditions
Equilibrium configurations in compression corresponding to the boundary conditions reported in Sect. 2.5 exhibit distinct height-radius profiles and thus different envelope surfaces, as shown in Fig. 6. However, all cases have an interesting common feature: solutions tend to overlap within a bulk region, while they depart from each other only at boundary layers, as one can appreciate from the insets in Fig. 6. Notably, the equilibrium solution for Case 3, corresponding to free radial displacement and radial rotations at terminal ends, results in a cylindrical deformation, i.e., a constant radius profile (see Sect. 2.3), which agrees with computational results in Ref. [11].

Curvatures and torsion
Curvatures and torsion profiles (Fig. 7) are symmetric with respect to the arc-length coordinate, consistently with the symmetry of the problems. They exhibit discontinuities at pin-joints, which correspond to concentrated moments exchanged between the helices. Figure 8 illustrates the balance of moments at pin-joints for the first right-handed helix, which is consistent with the absence of external couples and shows a symmetric pattern as well. All cases show a bulk region in the center of the assembly, where curvatures and torsion tend to be constant with similar values across cases, and two boundary layers, where they tend to differ (Fig. 7b) or to oscillate around the bulk value (Fig. 7a). Remarkably, Case 3 shows constant curvatures and torsion profiles, without boundary layers. Hence, in this case rods remain circular helices under compression, while the other cases approach the helical assembly of Case 3 in their bulk regions, and they   Figure 9 reports a comparison of the axial force versus compressive axial strain |ε| curves for the four cases treated, with ε = 0 and ε = −1 corresponding to the initial and a planar circle-like configurations, respectively. All cases show similar trends for the vertical force as a function of |ε|. Indeed, all force profiles are non-linear and monotonic, with higher forces for higher compressions and, for a given strain, they have similar values. As compression increases, the required external force grows, since it is balanced by a growing elastic restoring force.

Case of zero intrinsic curvatures and torsion with an example of zero-stiffness structure
In previous cases rods were modeled with linear constitutive equations, setting intrinsic curvatures and torsion equal to the initial ones. This choice implies that the initial equilibrium configuration is stress-free. A different scenario arises when κ s 1 = κ s 2 = Ω s = 0, that is, for naturally straight rods, since the initial configuration now requires forces and moments to be in equilibrium. Then, we can obtain the force-strain curve of a single circular helix by simplifying Eq. (48) aŝ Figure 10 compares the force profile for an assembly with κ s 1 = κ s 2 = Ω s = 0 and boundary conditions corresponding to Case 3 as given by numerical simulations, with that of Eq. (57), highlighting their perfect matching. Moreover, starting from this model we can envision a zero-stiffness structure. Indeed, according to Eq. (57), the vertical force is linear with respect to the normalized height h with a coefficient proportional to K 3 − 1, where K 3 is the ratio of torsional and bending stiffnesses. Therefore, by setting K 3 = 1 we obtain an assembly requiring a constant zero force for any compression ratio, as confirmed numerically (Fig. 10). Deformed configurations form an energetic plateau, where bending and torsional energies are interchangeable.

Comparison with Sadowsky's ribbon model
Using the constitutive Eqs. (19) and (20) and assuming a rectangular cross-section, Kirchhoff's rod model allows for a maximum stiffness for the force-strain relation (57) given by where β is a geometric parameter which increases with increasing aspect ratio of the cross-section, and that tends to 1/3 when the aspect ratio tends to infinity [25]. We therefore considered an alternative model for the elasticity of the rods, namely, the model for ribbon-like structures proposed by Sadowsky [26] and re-examined in Refs. [27][28][29], in order to investigate the impact of the choice of different constitutive models on the force-strain response. When the rods are modeled as narrow elastic ribbons, the energy density of each rod is described as where B = wt 3 12 with w the width of the ribbon, t its thickness, E Y its Young's modulus and ν its Poisson ratio. It is assumed that t w L, where L is the length of the rod, and that intrinsic intrinsic curvatures and torsion are equal to zero [27][28][29]. The energy density of Eq. (59) is not strictly convex (see Fig. 12), and this causes convergence problems in the numerical simulations. To circumvent this issue, we added a corrective term to the energy density and defined where ζ > 0 is a regularization parameter. The energy density of Eq. (61) is strictly convex for all ζ > 0, and Sadowski's energy (59) is recovered in the limit ζ → 0. The choice ζ = 0.2B guaranteed convergence in all the numerical simulations we conducted. Given the boundary conditions of Case 3, we simulated the extension of the assembly from h = h 0 = 0.01 L to h = L and we compared the results with those for a rod model assuming κ s = 0, Ω s = 0, ν = 0.3, and w t, which corresponds to K 3 ≈ 1.54, i.e., maximum stiffness according to Eq. (58). Both models adopt the same material and geometric parameters.
In the case of the ribbon model, to obtain the response of the assembly in the absence of the regularizing term we need a procedure to send ζ to zero. This is, however, easy when considering the boundary conditions of Case 3, because the equilibrium configurations that we obtain using the regularized energy, Eq. (61), are given by cylindrical deformations for any imposed height h. At a certain h, the same cylindrical deformation (modulo rigid body motions or change of chirality) is an equilibrium solution for any ζ > 0, and it is unique thanks to the convexity of the energy density. The axial forces experienced by rods are insensitive to ζ because the correction term in Eq. (61), when evaluated on circular helices, contributes with zero axial force (see Eq. (57), K 3 = 1). In these configurations, the (constant) internal moments are given by where we have written only the non-zero components. At each h, the same cylindrical deformation is an equilibrium solution also for the non-regularized (ζ = 0) energy density of Eq. (59), under the same axial force and internal moments given by M 1 (h, 0), M 3 (h, 0). This can be easily verified by checking that it satisfies the differential equilibrium equations, Eqs. (16) and (17). The axial force-axial strain profiles obtained using the two models are compared in Fig. 11. Since a ribbon is a shelllike structure with no stretches of the mid-surface allowed, while no such constraint is present in the rod model, it is not surprising that the ribbon model leads to a stiffer response of the assembly. What we find surprising, and do not fully understand, is the magnitude of the discrepancy in the prediction of the response using the two different models. It would be interesting to investigate this issue further, and to compare analytical predictions with experimental results, obtained for rods of rectangular cross-section with different aspect ratios w/t. One could then establish which model, if any, better captures the experimentally observed response, and in which regime of geometric or material parameters. Another interesting feature of the ribbon model is that it leads to a discontinuity in the tangent of the axial force-axial strain profile, i.e., in the axial stiffness of the assembly. This is caused by a discontinuity of the second derivative of the energy density, which occurs at |κ 1 | = |Ω|, as illustrated in Fig. 12 (ζ = 0). Indeed, in the region |κ 1 | < |Ω| the energy density of a rod becomes degenerate (independent of κ 1 ) and the energetic landscape is no longer strictly convex. A signature of this degeneracy is visible in the |κ 1 | − |Ω| trajectories traced by evaluating at the pin-joints of one rod the equilibrium solutions obtained by using Sadowski's energy density without regularizing terms, Eq. (59), in a numerical simulation of extension (h 0 = 0.01 L, see Fig. 12). As long as the solution remains in the strictly convex region |κ 1 | > |Ω|, we obtain circular helices (i.e., constant curvature and torsion) represented at each h by a single point (|κ 1 (h)|, |Ω(h)|) along the curve given by Eqs. (29) and (30), namely, κ 1 (ĥ) = 2 π n 1 −ĥ 2 ,Ω(ĥ) = 2 π nĥ .
(63) Fig. 11 Normalized axial force-axial strain profiles for a rod model with κ s = 0, Ω s = 0, K 3 ≈ 1.54, corresponding to maximum stiffness (orange), and for Sadowsky's ribbon model (blue). Insets compare the geometric configurations obtained in the two models at a certain level of axial force; the initial configuration is reported as a reference Fig. 12 Energetic landscape of Sadowsky's ribbon model (ζ = 0), over which we drew the |κ 1 |-|Ω| curves evaluated at the pin-joints of a rod, as traced by solutions of the equilibrium problem. The inset highlights their "fraying" in the non-convex region (|κ 1 | < |Ω|) As soon as the solution crosses the line |κ 1 | = |Ω|, the trajectories "fray" (that is, the single point at each h degenerates into a point cloud) indicating that the rod is exploring zero-energy perturbations, i.e., changes ofκ 1 at constantΩ (see the inset in Fig. 12). After a few steps, the simulation stops converging.

Conclusions and outlook
Focusing on a concrete example, we have shown that, when arranged in assemblies, helical shapes stabilize each other thanks to the interactions arising from mutual constraints. The pin-jointed helices of our example respond collectively to extension and compression, preserving the axial symmetry of the initial configuration. They can remain cylindrical, i.e., each rod can deform as a circular helix. This is the case of rods with ends free to displace, yet constrained to radial rotations. On the other hand, when ends are constrained from moving, the assembly forms edge layers, where rods behave like perturbed helices, and a bulk region, where circular helices are preserved. These circular helix configurations characterize the response of an assembly of infinite length [23]. Our results suggest that, boundary conditions incompatible with the geometry of a circular helix generate a perturbation in the form of boundary layers which act as transitions from imposed constraints to a bulk region where rods deform as circular helices. In future work, it would be interesting to study the dependence of the extension of the edge layers from the geometric and material parameters of the system, e.g., from the initial height of the assembly.
Moreover, we are interested in considering the use of different models for the prediction of the response of tubular helical structures, to overcome the limitations of the ones we have adopted so far. It would be interesting to experimentally validate the axial force-axial strain profiles in Fig. 11, in order to verify which model, and in which regime, best reproduces the physical behavior of the assembly, as it has been done for knots in Ref. [29]. Our observation on the ill-posedness of the Sadowsky's model, corrected as in Ref. [28], also raises the question of deriving less degenerate Sadowsky-type energies for narrow ribbons, in which higher order terms may provide a physically meaningful correction and lead to a less singular behavior.
At a more general level, the collective behavior we have discussed highlights the crucial role of mutual interactions among the rods in promoting the stability and persistence of helical shapes in a tubular assembly, hinting at why they are so recurrent in natural and artificial systems. The pellicle strips of Euglena shown in Fig. 1b, in which adjacent strips slide along their common edge, provide another example of stabilization of helical shapes due to mutual interactions in a tubular assembly [9,30]. It would be interesting to see whether this argument can be extended further to other, less specific types of interactions among the fibers of the assembly, such as in the microtubule meshwork of Lacrymaria shown in Fig. 1a.
Funding Open access funding provided by Scuola Superiore Sant'Anna within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.