Forward and inverse problems in the mechanics of soft filaments

Soft slender structures are ubiquitous in natural and artificial systems, in active and passive settings and across scales, from polymers and flagella, to snakes and space tethers. In this paper, we demonstrate the use of a simple and practical numerical implementation based on the Cosserat rod model to simulate the dynamics of filaments that can bend, twist, stretch and shear while interacting with complex environments via muscular activity, surface contact, friction and hydrodynamics. We validate our simulations by solving a number of forward problems involving the mechanics of passive filaments and comparing them with known analytical results, and extend them to study instabilities in stretched and twisted filaments that form solenoidal and plectonemic structures. We then study active filaments such as snakes and other slender organisms by solving inverse problems to identify optimal gaits for limbless locomotion on solid surfaces and in bulk liquids.

we demonstrate the capabilities and the robustness of our solver by embedding it in an inverse design cycle for the identification of optimal terrestrial and aquatic limbless locomotion strategies.
The paper is structured as follows. In §2, we review and introduce the mathematical foundations of the model. In §3, we present the corresponding discrete scheme and validate our implementation against a battery of benchmark problems. In §4, we detail the physical and biological enhancements to the original model, and finally in §5, we showcase the potential of our solver via the study of solenoids and plectonemes as well as limbless biolocomotion. Mathematical derivations and additional validation test cases are presented in the appendix.

Governing equations
We consider filaments as slender cylindrical structures deforming in three dimensions with a characteristic length L which is assumed to be much larger than the radius (L r) at any cross section. Then the filament can be geometrically reduced to a one-dimensional representation, and its dynamical behaviour may be approximated by averaging all balance laws at every cross section [5]. We start with a description of the commonly used Kirchhoff-Love theory that accounts for bend and twist at every cross section but ignores stretch and shear, before moving on to the Cosserat theory that also accounts for these additional degrees of freedom.

Kirchhoff-Love theory for inextensible, unshearable rods
As illustrated in figure 1, a filament in the Cosserat rod theory can be described by a centre line r : (s ∈ [0, L] ∈ R, t ∈ R + ) → R 3 and an oriented frame of reference Q : (s ∈ [0, L] ∈ R, t ∈ R + ) → SO(3) equivalent to the orthonormal triad of unit vectors Q = {d 1 , d 2 , d 3 }. Here, s is the centre line arc-length coordinate in its current configuration and t is the time.
Denoting by x any generic vector represented in the Eulerian frame and x L as the body-convected (Lagrangian) frame of reference allows us to write laboratory: x =x 1 i +x 2 j +x 3 k (2.1) and body-convected: where equation (2.1) expresses x in the laboratory canonical basis {i, j, k}, while equation (2.2) expresses the same vector in the body-convected director basis {d 1 , d 2 , d 3 }. Then, the matrix Q transforms any vector x from the laboratory to the body-convected representation via x L = Qx and, conversely, x = Q −1 x L = Q T x L , because Q T Q = QQ T = 1 1 1. In general, we need to distinguish between the arc-length coordinate s that corresponds to the current filament configuration and the arc-length coordinateŝ associated with the reference configuration of the filament due to stretching (figure 1-throughout we will use hatted quantities to denote the reference configuration). We will first start by presenting the equations of motion under the assumption of inextensibility (i.e. s =ŝ, Kirchhoff model), before generalizing them to the stretchable case in the subsequent sections. Denoting the rod angular velocity as ω = vec[(∂Q/∂t) T Q] and the generalized curvature as κ = vec[(∂Q/∂s) T Q], where vec[A] denotes the 3-vector associated with the skew-symmetric matrix A, the following transport identities hold Using the above equations (full derivation in the appendix), we can express the advection of the rod positions and local frames, as well as the linear and angular momentum balance in a convenient Eulerian-Lagrangian form  Figure 1. The Cosserat rod model. A filament deforming in the three-dimensional space is represented by a centre-line coordinate r and a material frame characterized by the orthonormal triad {d 1 , d 2 , d 3 }. The corresponding orthogonal rotation matrix Q with row entries d 1 , d 2 , d 3 transforms a vector x from the laboratory canonical basis {i, j, k} to the material frame of reference {d 1 , d 2 , d 3 } so that x L = Qx and vice versa x = Q T x L . If extension or compression is allowed, the current filament configuration arc-length s may no longer coincide with the rest reference arc-lengthŝ. This is captured via the scalar dilatation field e = ds/dŝ. Moreover, to account for shear we allow the triad {d 1 , d 2 , d 3 } to detach from the unit tangent vector t so that d 3 = t (we recall that the condition d 3 = t and e = 1 correspond to the Kirchhoff constraint for unshearable and inextensible rods, and implies that σ = et − d 3 = 0). The dynamics of the centre line and the material frame is determined at each cross section by the internal force and torque resultants, n and τ . External loads are represented via the force f and couple c line densities. Table 1. Constitutive laws. The generalized curvature κ L is associated with the bending κ 1 , κ 2 about the principal directions (d 1 , d 2 ) and the twist κ 3 about the longitudinal one (d 3 ), while σ L = Q(et − d 3 ) is associated with the shears σ 1 , σ 2 along the principal directions (d 1 , d 2 ) and the axial extensional or compression σ 3 along the longitudinal one (d 3 ). The material properties of the rod are captured through the Young's (E) and shear (G) moduli, while its geometric properties are accounted for via the cross-sectional area A, the second moment of inertia I and the constant α c = 4 3 for circular cross sections [58]. The diagonal entries of the bending/twist B ∈ R 3×3 and shear/stretch S ∈ R 3×3 matrices are, respectively, (B 1  where ρ is the constant material density, A is the cross-sectional area, v is the velocity, n L and τ L are, respectively, the internal force and couple resultants, f and c are external body force and torque line densities, and the tensor I is the second area moment of inertia (throughout this study we assume circular cross sections; see the appendix).
To close the above system of equations (2.4)-(2.7) and determine the dynamics of the rod, it is necessary to specify the form of the internal forces and torques generated in response to bend and twist, corresponding to the 3 d.f. at every cross section. The strains are defined as the relative local deformations of the rod with respect to its natural strain-free reference configuration. Bending and twisting strains are associated with the spatial derivatives of the material frame director field {d 1 , d 2 , d 3 } and are characterized by the generalized curvature. Specifically, the components of the curvature projected along the directors (κ L = κ 1 d 1 + κ 2 d 2 + κ 3 d 3 ) coincide with bending (κ 1 , κ 2 ) and twist (κ 3 ) strains in the material frame (table 1) We assume a perfectly elastic material so that the stress-strain relations are linear. Integration of the torque densities over the cross-sectional area A yields the bending and twist rigidities (table 1), so that the resultant torque-curvature relations can be generically expressed in vectorial notation as is the bend/twist stiffness matrix, with B 1 the flexural rigidity about d 1 , B 2 the flexural rigidity about d 2 and B 3 the twist rigidity about d 3 . Here, the vector κ o L characterizes the intrinsic curvatures of a filament that in its stress-free state is not straight. We wish to emphasize here that the constitutive laws are most simply expressed in a local Lagrangian form; hence the use of κ L and not κ.
The Kirchhoff rod is defined by the additional assumption that there is no axial extension or compression or shear strain. Then the arc-length s coincides withŝ at all times, and the tangent to the centre line is also normal to the cross section, so that t = d 3 [5]. This implies that n L serves as a Lagrange multiplier, and that the torque-curvature relations of equation (2.8) are linear.
This completes the formulation of the equations of motion for the Kirchhoff rod, and when combined with boundary conditions suffices to have a well-posed initial boundary value problem. For the general stretchable and shearable case, all geometric quantities (A, I, κ L , etc.) must be rescaled appropriately, as addressed in the following sections.

Cosserat theory of stretchable and shearable filaments
In the general case of soft filaments, at every cross section we also wish to capture transverse shear and axial strains in addition to bending and twisting. As we wish to account for all six deformation modes associated with the 6 d.f. at each cross section along the rod, we must augment the Kirchhoff description of the previous section and add three more constitutive laws to define the local stress resultants n L . In fact, they are no longer defined as Lagrange multipliers that enforce the condition of inextensibility and unshearability, i.e. we now must allow t = d 3 . We note that this not only enables the model of equations (2.4)-(2.7) to capture a richer set of physical phenomena, but also significantly simplifies its numerical treatment (Lagrange multipliers no longer needed), rendering this model both flexible and relatively straightforward to implement.
The shear and axial strains are associated with the deviations between the unit vector perpendicular to the cross section and the tangent to the centre line, and thus may be expressed in terms of the derivatives of the centre line coordinate r. In the material frame of reference, we characterize these strains by the vector σ (figure 1) which then takes the form Here, the scalar field e(ŝ, t) = ds/dŝ expresses the local stretching or compression ratio (figure 1) relative to the rest reference configuration (ŝ) and t is the unit tangent vector. Whenever the filament undergoes axial stretching or compression the corresponding infinitesimal elements deform and all related geometric quantities are affected. By assuming that the material is incompressible and that the cross sections retain their circular shapes at all times, we can conveniently express the governing equations with respect to the rest reference configuration of the filament (denoted by a hat) in terms of the local dilatation e(ŝ, t). Then, the following relations hold: (2.10) As with the Kirchhoff rod, assuming a linear material constitutive law implies linear stress-strain relations. Integration of the stress and couple densities over the cross-sectional area A yields both the rigidities associated with axial extension and shear (table 1), so that the resultant load-strain relations can be generically expressed in vectorial notation as the vector σ o L corresponds to the intrinsic shear and stretch, and must be accounted for in the case of stress-free, non-trivial shapes. Although the intrinsic strains κ o L , σ o L are implemented in our solver to account for pre-strained configurations, to simplify the notation in the remaining text, we will assume that the filament is intrinsically straight in a stress-free state, so that σ o L = κ o L = 0. The rigidities associated with bending, twisting, stretching and shearing are specified in table 1, and can be expressed as the product of a material component, represented by the Young's (E) and shear (G) moduli, and a geometric component represented by A, I and the constant α c = 4 3 for circular cross sections [58]. We note that the rigidity matrices B and S are assumed to be diagonal throughout this study, although off-diagonal entries can be easily accommodated to model anisotropic materials such as composite elements. In general, this mathematical formulation can be extended to tackle a richer set of physical problems including viscous threads [15,16], magnetic filaments [59], etc., by simply modifying the entries of B and S and introducing time-dependent constitutive laws wherein τ L (κ L , ∂ t κ L ) and n L (σ L , ∂ t σ L ), as for example in [15]. We also emphasize here that in the case of stretchable rods, A and I are no longer constant, rendering the load-strain relations nonlinear (equations (2.8) and (2.11)), even though the stress-strain relations remain linear. This is consistent with the modelling of hyperelastic materials such as rubbers, silicones and biological tissues and therefore in line with targeted soft robotic applications. Indeed, the combination of linear stress-strain constitutive models with the geometrical rescaling by e leads to a reasonable approximation of Neo-Hookean [60] and Mooney-Rivlin models [61,62] (especially tailored to hyperelastic solids) over a compression/extension range up to 30%. See the appendix for a quantitative comparison in the axial stretch case.
Having generalized the constitutive relations to account for filament stretchiness and shearability, we now generalize the equations of motion for this case. Multiplying both sides of equations (2.6) and (2.7) by ds and substituting the above identities together with the constitutive laws of equation (2.11) into equations (2.4)-(2.7) yields the final system where dm = ρÂ dŝ = ρA ds is the infinitesimal mass element, and dĴ = ρÎ dŝ is the infinitesimal mass second moment of inertia. We note that the left-hand side and the unsteady dilatation term of equation (2.15) arise from the expansion of the original rescaled angular momentum dĴ · ∂ t (ω L /e) via the chain rule. We also note that the external force and couple are defined as F = ef dŝ and C L = ec L dŝ (with f and c L the force and torque line densities, respectively) so as to account for the dependence on e. Combined with some initial and boundary conditions, equations (2.12)-(2.15) express the dynamics and kinematics of the Cosserat rod with respect to its initial rest configuration, in a form suitable to be discretized as described in §3.

Numerical method
To derive the numerical method for the time evolution of a filament in analogy with the continuum model of §2, we first recall a few useful definitions for effectively implementing rotations. We then present the spatially discretized model of the rod, and the time discretization approach employed to evolve the governing equations.  Figure 2. Discretization model. A discrete filament is represented through a set of vertices r(t) i=1,...,N+1 and a set of material frames ..,N . Two consecutive vertices define an edge of length i along the tangent unit vector t i . The dilatation is defined as e i = i /ˆ i , whereˆ i is the edge rest length. The vector σ i = e i t i − d i 3 represents the discrete shear and axial strains. The mass m r of the filament is discretized in pointwise concentrated masses m i=1,...,n+1 at the locations r i for the purpose of advecting the vertices in time. For the evolution of Q i in time, we consider instead the mass second moment of inertiaĴ i=1,...,n associated with the cylindrical elements depicted in blue.

Rotations
Bending and twisting deformations of a filament involve rotations of its material frame Q in space and time. To numerically simulate the rod, it is critical to represent and efficiently compute these nonlinear geometric transformations quickly and accurately. A convenient way to express rotations in space or time is the matrix exponential [63][64][65]. Assuming that the matrix R denotes the rotation by the angle θ about the unit vector axis u, this rotation can be expressed through the exponential matrix R = e θu : R 3 → R 3×3 , and efficiently computed via the Rodrigues formula [66] e θu = 1 1 1 + sin θU + (1 − cos θ)UU. (3.1) Here, U ∈ R 3×3 represents the skew-symmetric matrix associated with the unit vector u where the operator [·] × : R 3 → R 3×3 allows us to transform a vector into the corresponding skewsymmetric matrix, and vice versa [·] −1 × : R 3×3 → R 3 . Conversely, given a rotation matrix R, the corresponding rotation vector can be directly computed via the matrix logarithm operator log(·) : It is important to note that the rotation axis u is expressed in the material frame of reference associated with the matrix R (or Q). With these tools in hand, we now proceed to outline our numerical scheme.

Spatial discretization
Drawing from previous studies of unshearable and inextensible rods [13,14,67], we capture the deformation of a filament in three-dimensional space via the time evolution of a discrete set of vertices r i (t) ∈ R 3 , i ∈ [1, n + 1] and a discrete set of material frames Q i (t) ∈ R 3×3 , i ∈ [1, n], as illustrated in figure 2.
Each vertex is associated with the following discrete quantities: where v i is the velocity , m i is a pointwise concentrated mass and F i is the external force given in equation (2.14). Each material frame is associated with an edge i connecting two consecutive vertices, and with the related discrete quantities where i = | i |,ˆ i = |ˆ i |, e i = i /ˆ i are the edge current length, reference length and dilatation factor, t i is the discrete tangent vector, σ i L is the discrete shear/axial strain vector, ω i L is the discrete angular velocity,Â i ,Ĵ i ,B i ,Ŝ i are the edge reference cross-section area, mass second moment of inertia, bend/twist matrix and shear/stretch matrix, respectively, and finally C i L is the external couple given in equation (2.15).
Whereas in the continuum setting ( §2) all quantities are defined pointwise, in a discrete setting some quantities, and in particular κ L , are naturally expressed in an integrated form over the domain D along the filament [14,68]. Any integrated quantity divided by the corresponding integration domain length D = |D| is equivalent to its pointwise average. Therefore, following the approach of Bergou & Grinspun [14,68], the domain D becomes the Voronoi region D i of length which is defined only for the interior vertices r ,...,n−1 . Each interior vertex is then also associated with the following discrete quantities: whereD i is the Voronoi domain length at rest and E i is Voronoi region dilatation factor. Recalling that the generalized curvature expresses a rotation per unit length about its axis, the quantityD iκ i L naturally expresses the rotation that transforms a material frame Q i to its neighbour Q i+1 over the segment sizê D i along the rod. Therefore, the relation eD iκ i L Q i = Q i+1 holds, so thatκ i L = log(Q i+1 Q T i )/D i . Finally, we introduce the bend/twist stiffness matrixB i consistent with the Voronoi representation.
Then, we may discretize the governing equations (2.12)-(2.15) so that they read where h : is the averaging operator (trapezoidal quadrature rule) to transform integrated quantities over the domain D to their point-wise counterparts [14,69]. We note that h and A h operate on a set of N vectors and returns N + 1 vectors, consistent with equations (3.6)-(3.9) (see the appendix for further details).

Time discretization
The derivation above leads to a system that in general is not Hamiltonian, as this depends on the nature of the external loads acting on the filament as well as on the choice of constitutive laws. Nonetheless, without external loads and under the assumptions of linear stress-strain relations, this derivation amounts to a geometric rescaling through e of the classic Cosserat rod model (here directly discretized via standard finite difference and trapezoidal quadrature rule, as outlined in [14,18,69,70]), which is a Hamiltonian system [18,71] with quadratic energy functionals [70]: Therefore, by construction, in the limit of e → 1 (no axial elongation) the equivalence with the Cosserat model holds, and for consistency we opt for an energy-preserving time integrator, and in particular a symplectic, second-order Verlet scheme. We note that, despite the failure of Verlet schemes to integrate rotational equations of motion when represented by quaternions, in our case their use is justified as rotations are represented instead by Euler angles [72].
The second-order position Verlet time integrator is structured in three blocks: a first half-step updates the linear and angular positions, followed by the evaluation of local linear and angular accelerations, and finally a second half-step updates the linear and angular positions again. Therefore, it entails only one right-hand side evaluation of equations (3.8) and (3.9), the most computationally expensive operation (see the appendix for details).
This algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity: by foregoing an implicit integration scheme we can incorporate a number of additional physical effects and soft constraints, even though this may come at the expense of computational efficiency. Indeed, for large Young's or shear moduli or for very thin rods, the system of equations (2.4)-(2.7) might become stiff, so that small timesteps must be employed to ensure stability. Although we have not derived a rigorous CFL-like condition, throughout this work we employed the empirical relation dt = a d , with a ∼ 10 −2 s/m, and found it reliable in preventing numerical instabilities. Despite the potential stiffness of this model, as the computational cost per timestep scales linearly with the number of discretization elements n (the model is one-dimensional), we could carry out all our computation on a laptop (details and representative timings are summarized in the appendix, figure 20). It is worth noting that because of the condition dt = a d , a smaller number of elements implies larger timesteps (d = L/n). Therefore, the time-to-solution may scale between linearly and quadratically, depending on how the timestep is set.

Validation
We first validate our proposed methodology against a number of benchmark problems with analytic solutions and examine the convergence properties of our approach. Three case studies serve to characterize the competition between bending and twisting effects in the context of helical buckling, dynamic stretching of a loaded rod under gravity, and the competition between shearing and bending in the context of a Timoshenko beam. Further validations reported in the appendix include Euler and Mitchell buckling due to compression or twist, and stretching and twisting vibrations.

Helical buckling instability
We validate our discrete derivative operators beyond the onset of instability (see Euler and Mitchell buckling tests in the appendix) for a long straight, isotropic, inextensible, and unshearable rod undergoing bending and twisting. The filament is characterized by the length L and by the bending and twist stiffnesses α and β, respectively. The clamped ends of the rod are pulled together in the axial direction k with a slack D/2 and simultaneously twisted by the angle Φ/2, as illustrated in figure 3a. Under these conditions, the filament buckles into a localized helical shape (figure 3e).
The nonlinear equilibrium configuration r eq of the rod can be analytically determined [8,[73][74][75] in terms of the total applied slack D and twist Φ. We denote the magnitude of the twisting torque and tension acting on both ends and projected on k by M h and T h , respectively. Their normalized counterparts m h = M h L/(2πα) and t h = T h L 2 /(4π 2 α) can be computed via the 'semi-finite' correction approach [74] by solving the system . (e) Equilibrium rod configuration r n eq numerically obtained given the discretization n = 800, and assuming dissipation. For all studies, unless specified otherwise, we used the following settings: length L = 100 m, twist Φ = 27 · 2π , slack D = 3 m, linear mass density ρ = 1 kg m −1 , bending stiffness α = 1.345 Nm 2 , twisting stiffness β = 0.789 Nm 2 , shear/stretch matrix S = 10 5 · 1 1 1 N, bend/twist matrix B = diag(α, α, β) Nm 2 , dissipation constant γ = 10 −2 kg (ms) −1 , radius r = 0.35 m, twisting time T twist = 500 s, relaxation time T relax = 10 4 s.
Then, the analytical form of r eq can be expressed [75] as wheres = s/L − 0.5 is the normalized arc-length −0.5 ≤s ≤ 0.5. Here, we make use of equation (3.10) to investigate the convergence properties of our solver in the limit of refinement. To compare analytical and numerical solutions, a metric invariant to rotations about k is necessary. Following Bergou et al. [14], we rely on the definition of the envelope ϕ ϕ = cos θ − cos θ max 1 − cos θ max , θ = arccos(t · k) (3.11) where θ is the angular deviation of the tangent t from the axial direction k, and θ max is the corresponding maximum value along the filament. The envelope ϕ relative to the analytical solution of equation (3.10), and ϕ n relative to a numerical model of n discretization elements can be estimated via finite differences. This allows us to determine the convergence order of the solver by means of the norms L 1 ( ), L 2 ( ) and L ∞ ( ) of the error = ϕ − ϕ n . We simulate the problem illustrated in figure 3 at different space-time resolutions. The straight rod originally at rest is twisted and compressed at a constant rate during the period T twist . Subsequently, the ends of the rod are held in their final configurations for the period T relax to allow the internal energy to dissipate (according to the model of §4.1) until the steady state is reached. Simulations are carried out progressively refining the spatial discretization δl = L/n by varying n = 100-3200 and the time discretization δt is kept proportional to δl, as reported in figure 3.
As can be seen in figure 3b,c,e the numerical solutions converge to the analytical one with second order in time and space, consistent with our spatial and temporal discretization schemes. Moreover, to validate the energy-conserving properties of the solver in the limit of e → 1, we turn off the internal dissipation (figure 3d) and observe that the total energy of the filament E F is constant after T twist and matches its

Vertical oscillations under gravity
We consider a system in which a rod hanging from one end and subject to gravity g oscillates due to a mass m p suspended at the other end, and due to its own mass m r , as depicted in figure 4a,d. This system is analogous to a mass-spring oscillator. The static solution is then obtained by integrating the infinitesimal elongations along the spring due to the local load [76], yielding the total equilibrium extension where k is the spring constant, ξ = 2 is a constant factor and m eq = m p + m r /ξ is the equivalent mass. Thus, the final equilibrium length of the rod reads L =L + L * , withL being the rest unstretched length. The dynamic solution is instead characterized by oscillations of period T * and by the time-varying length L(t) of the spring In this case, unlike the static solution, the factor ξ depends on the ratio m r /m p . In fact it can be shown [76] that ξ 3 for m r /m p → 0, and ξ π 2 /4 for m r /m p → ∞.
The analytical results rely on the assumption of k being constant in space and time, given a fixed ratio m r /m p . However, this condition is not met here because k(s, t) = EA(s, t) is a function of space and time, due to dilatation and mass conservation. Nevertheless, as the Young's modulus E → ∞, that is as a soft filament becomes stiff, the constant k → EÂ and our rod model must recover the behaviour of the massspring oscillator. Indeed, figure 4b,c,e,f shows how the proposed numerical method converges to the analytical oscillation period T * and normalized longitudinal displacement (L −L)/ L * as E increases.

Cantilever beam
We now consider the effect of bend and shear simultaneously by validating our numerical methods against the Timoshenko cantilever of figure 5a. Timoshenko's model accounts for bending elasticity, rotary inertia and shear deformations, building on classical beam theories by Rayleigh (bending elasticity and rotary inertia) and Euler-Bernoulli (bending elasticity only). The model captures the behaviour of short or composite beams in which shear deformations effectively lower the stiffness of the rod [58,77].
fi rs t o rd e r We consider a beam clamped at one endŝ = 0 and subject to the downward force F at the free end s =L, as illustrated in figure 5a. The static solution for the displacement y along the vertical direction i of the rod can then be analytically expressed as whereL is the length of the rod,Â is the constant cross-sectional area,Î 1 is the area second moment of inertia about the axis j = k × i, E and G are, respectively, the Young's and shear moduli, and α c = 4 3 is the Timoshenko shear factor for circular sections and accounts for the fact that the shear stress varies over the section [58]. Furthermore, the Timoshenko (as well as Rayleigh and Euler-Bernoulli) theory relies on the assumption of small deflections, so that the horizontal coordinate x along the direction k can be approximated by the arc-lengthŝ (see figure 5a and the appendix for further details and derivation), hence the use ofŝ in the above equation.
If the shear modulus G approaches infinity or if the ratio EÎ 1 /(α cL 2Â G) 1, then the Timoshenko model in the static case reduces to the Euler-Bernoulli approximation, yielding as the shear term of equation (3.14) becomes negligible.
We compare our numerical model with these results by carrying out simulations of the cantilever beam of figure 5a in the time-space limit of refinement. As can be noticed in figure 5b the discrete solution recovers the Timoshenko one. Therefore, the solver correctly captures the role of shear that reduces the effective stiffness relative to the Euler-Bernoulli solution. Moreover, our approach is shown to converge to the analytical solution in all the norms L ∞ ( ), L 1 ( ), L 2 ( ) of the error = y − y n , where y n is the vertical displacement numerically obtained in the refinement limit.
We note that the norms L ∞ ( ) and L 1 ( ) exhibit first-order convergence, while L 2 ( ) decays with a slope between first and second order. We attribute these results to the fact that while the Timoshenko solution does not consider axial extension or tension, it does rely on the assumption of small deflections (ŝ = x), therefore effectively producing a dilatation of the rod. On the contrary, our solver does not assume small deflections and does not neglect axial extension, because the third entry of the matrix B has the finite value EÂ (see figure 5 for details). This discrepancy is here empirically observed to decrease the convergence order.
These studies, together with the ones reported in the appendix, complete the validation of our numerical implementation and demonstrate the accuracy of our solver in simulating soft filaments in simple settings.

Including interactions and activity: solid and liquid friction, contact and muscular effects
Motivated by advancements in the field of soft robotics [41,44,45], we wish to develop a robust and accurate framework for the characterization and computational design of soft slender structures interacting with complex environments. To this end, we expand the range of applications of our formalism by including additional physical effects, from viscous hydrodynamic forces in the slenderbody limit and surface solid friction to self-contact and active muscular activity. As a general strategy, all new external physical interactions are accounted for by lumping their contributions into the external forces and couples F and C L on the right-hand side of the linear and angular momentum balance equations (2.14) and ( Figure 6. Muscular activity. Example of muscular activity amplitude profile (solid black line) described by cubic B-spline through N m = 8 control points (Ŝ i , β i ) with i = 0, . . . , N m − 1. The control points are located along the filament at the positionsŜ i , and are associated with the amplitude values β i . The first and last control points are fixed so that (Ŝ 0 , β 0 ) = (0, 0) and (Ŝ N m −1 , β N m −1 ) = (L, 0), therefore assuming the ends of the deforming body to be free.
captured by adding their contributions directly to the internal force n L and torque τ L resultants before integrating equations (2.14) and (2.15).

Dissipation
Real materials are subject to internal friction and viscoelastic losses, which can be modelled by modifying the constitutive relations so that the internal torques τ L (κ L ) and forces n L (σ L ) of equation (2.11) become functions of both strain and rate of strain, i.e. τ L (κ L , ∂ t κ L ) and n L (σ L , ∂ t σ L ). Keeping track of the strain rates increases computational costs and the memory footprint of the solver. However, for the purpose of purely dissipating energy, a simple alternative option is to employ Rayleigh potentials [16,78]. In this case, viscous forces f v and torques c v L per unit length are directly computed as linear functions of linear and angular velocities through the constant translational γ t and rotational γ r internal friction coefficients, so that and This approach does not model the physical nature of viscoelastic phenomena, although it does dissipate energy, effectively mimicking overall material friction effects. In the context of our numerical investigations, we did not observe any appreciable difference between the two outlined methods, so that, for the sake of simplicity and computational efficiency, we opted for the second one. Throughout the remainder of the text we will then employ equations (4.1) and (4.2) with a single dissipation constant γ , therefore assuming γ t = γ r .

Muscular activity
To study limbless biolocomotion on solid substrates and in fluids, we allow our soft filaments to be active, by generating internal forces and torques corresponding to coordinated muscular activity driven, for example, by a central pattern generator [79,80].
Following the approach detailed in [53,56], we express the muscular activity magnitude A m as a travelling wave propagating head to tail along the filament where φ m is the phase, t is time and T m and λ m are, respectively, the activation period and wavelength. The amplitude of the travelling wave is represented by the cubic B-spline β(ŝ) characterized by N m control points (Ŝ i , β i ) with i = 0, . . . , N m − 1, as illustrated in figure 6. The first and last control points are fixed so that (Ŝ 0 , β 0 ) = (0, 0) and (Ŝ N m −1 , β N m −1 ) = (L, 0), therefore assuming the ends of the deforming body to be free. One of the main advantages of the proposed parametrization is that it encompasses a large variety of patterns with a relatively small number of parameters [56]. A given activation mode can be achieved by multiplying the scalar amplitude A m with the appropriate director. For example, if we wish to study earthworm-like locomotion, we may employ a wave of longitudinal dilatation and compression forces, so that assuming d 2 and d 3 to be coplanar to the ground. These two contributions are directly added to the internal force n L and torque τ L resultants.
In the most general case, all deformation modes can be excited by enabling force and torque muscular activity along all directors d 1 , d 2 and d 3 , providing great flexibility in terms of possible gaits.

Self-contact
To prevent the filament from passing through itself, we need to account for self-contact. As a general strategy, we avoid enforcing the presence of boundaries via Lagrangian constraints as their formulation may be cumbersome [81], impairing the modularity of the numerical solver. We instead resort to calculating forces and torques directly and replacing hard constraints with 'soft' displacement-force relations.
Our self-contact model introduces additional forces F sc acting between the discrete elements in contact. To determine whether any two cylindrical elements are in contact, we calculate the minimum distance d ij min between edges i, j by parametrizing their centre lines If d ij min is smaller than the sum of the radii of the two cylinders, then they are considered to be in contact and penalty forces are applied to each element as a function of the scalar overlap ij = (r i + r j − d ij min ), where r i and r j are the radii of edges i and j, respectively. If ij is smaller than zero, then the two edges are not in contact and no penalty is applied. Denoting as d ij min the unit vector pointing from the closest point on edge i to the closest point on edge j, the self-contact repulsion force is given by where H( ij ) denotes the Heaviside function and ensures that a repulsion force is produced only in case of contact ( ij ≥ 0). The first term within the square brackets expresses the linear response to the interpenetration distance as modulated by the stiffness k sc , while the second damping term models contact dissipation and is proportional to the coefficient γ sc and the interpenetration velocity v i − v j .

Contact with solid boundaries
To investigate scenarios in which filaments interact with the surrounding environment, we must also account for solid boundaries. By implementing the same approach outlined in the previous section, obstacles and surfaces are modelled as soft boundaries allowing for interpenetration with the elements of the rod (figure 7). The wall response F w ⊥ balances the sum of all forces F ⊥ that push the rod against the wall, and is complemented by the other two components which help prevent possible interpenetration due to numerics. The interpenetration distance triggers a normal elastic response proportional to the stiffness of the wall, while a dissipative term related to the normal velocity component of the filament with respect to the substrate accounts for a damping force, so that the overall wall response reads where H( ) denotes the Heaviside function and ensures that a wall force is produced only in case of contact ( ≥ 0). Here, u w ⊥ is the boundary outward normal (evaluated at the contact point, that is the contact location for which the normal passes through the centre of mass of the element), and k w and γ w are, respectively, the wall stiffness and dissipation coefficients.

Isotropic and anisotropic surface friction
Solid boundaries also affect the dynamics of the filament through surface friction, a complex physical phenomenon in which a range of factors are involved, from roughness and plasticity of the surfaces in contact to the kinematic initial conditions and geometric set-up. Here, we adopt the Amonton-Coulomb model, the simplest of friction models.
This model relates the normal force pushing a body onto a substrate to the friction force through the kinetic μ k and static μ s friction coefficients, depending on whether the contact surfaces are in relative motion or not.
Despite the simplicity of the model, its formulation and implementation may not necessarily be straightforward, especially in the case of rolling motions. Given the cylindrical geometry of our filaments,  The surface normal u w ⊥ determines the direction of the wall's response F w ⊥ to contact. We note that F w ⊥ balances the sum of all forces F ⊥ that push the rod against the wall, and is complemented by the other two components, which allow it to amend to possible interpenetration due to numerics. These components are an elastic one (k w ) and a dissipative one (γ w v · u w ⊥ ), where k w and γ w are, respectively, the wall stiffness and dissipation coefficients. ⊥ and a longitudinal one in the direction u w × = u w ⊥ × u w . We note that in general u w × = t. The notation x ⊥ , x , x × denotes the projection of the vector x in the directions u w ⊥ , u w , u w × . (b,c) Kinematic and dynamic quantities at play at any cross section in the case of (b) rolling and slipping and (c) pure rolling motion. Red arrows correspond to forces and torques, green arrows correspond to velocities, and black arrows correspond to geometric quantities. the effect of surface friction can be decomposed into a longitudinal component associated with purely translational displacements, and a lateral component associated with both translational and rotational motions (figure 8). We use the notation x ⊥ , x , x × to denote the projection of the vector x in the directions u w ⊥ , u w , u w × , as illustrated in figure 8. The longitudinal friction force F long is opposite to either the resultant of all forces F × acting on an element (static case) or to the translational velocity v × (kinetic case) along the direction u w × ( figure 8). The Amonton-Coulomb model then reads where v → 0 is the absolute velocity threshold value employed to distinguish between static (|v × | ≤ v ) and kinetic (|v × | > v ) cases. We define v in a limit form to accommodate the fact that inequalities are numerically evaluated up to a small threshold value. The static friction force is always equal and opposite to F × up to a maximum value proportional to the normal force |F ⊥ | through the coefficient μ s . The kinetic friction force is instead opposite to the translational velocity v × , but does not depend on its actual magnitude and is proportional to |F ⊥ | via μ k . In general μ s > μ k , so that it is harder to set a body into motion from rest than to drag it. The lateral displacement of a filament in the direction u w = u w × × u w ⊥ is associated with both translational (v ) and rotational (ω × = ω × u w × ) motions, as illustrated in figure 8b,c. In this case, the distinction between static and kinetic friction does not depend on v , but on the relative velocity (also referred to as slip velocity) between the rod and the substrate where v cont is the local velocity of the filament at the contact point with the substrate, due to the axial component of the angular velocity ω × .
In the static or no-slip scenario (v slip = 0), the linear momentum balance in the direction u w , and the angular momentum balance about the axis u w × express a kinematic constraint between the linear acceleration au w and angular acceleration ω × = (u w ⊥ × au w )/r, so that (F + F roll )u w = dm · au w (4.10) and where F = F u w and T × = T × u w × are the forces and torques acting on the local element, and F roll = F roll u w is the rolling friction force at the substrate-filament interface necessary to meet the no-slip condition. By recalling that a disk mass second moment of inertia about u w × is J = r 2 dm/2, the above system can be solved for the unknown a and F roll , yielding Therefore, the lateral friction force F lat and the associated torque C lat L can be finally expressed as where μ r s and μ r k are, respectively, the rolling static and kinetic friction coefficients. So far we have considered isotropic friction by assuming that the coefficients μ s and μ k are constant and independent of the direction of the total acting forces (static case) or relative velocities (kinetic case). Nevertheless, frictional forces may be highly anisotropic. For example, the anisotropy caused by the presence of scales on the body of a snake crucially affects gaits and performance [82,83].
The Amonton-Coulomb model can be readily extended to account for anisotropic effects by simply assuming the friction coefficients μ s and μ k to be functions of a given reference direction. The nature of these functions depends on the specific physical problems under investigation. An example of this approach is illustrated in §5.2 in the context of limbless locomotion. Isotropic and anisotropic friction validation benchmarks are presented in the appendix.

Hydrodynamics
We also extend our computational framework to address flow-structure interaction problems. In particular, we consider the case in which viscous forces dominate over inertial effects, i.e. we consider systems in which the Reynolds number Re = ρ f UL/μ 1 where ρ f and μ are the density and dynamic viscosity of the fluid, respectively, and U is the characteristic velocity of the rod. Under these conditions, the drag forces exerted by the fluid on our filaments can be determined analytically within the context of slender-body theory [84,85] (for more advanced and accurate viscous flow models we refer the reader to [31][32][33][34][35][36]). At leading order resistive force line densities scale linearly with the local rod velocities v according to We note that the matrix (I − 1 2 t T t) introduces an anisotropic effect for which (4.14) where f H = (f H · t)t and f H ⊥ = f H − f H are, respectively, tangential and orthogonal viscous drag components. The coupling of liquid environment, filament mechanics and muscular activity provides a flexible platform to characterize biological locomotion at the microscopic scale (bacteria, protozoa, algae, etc.) and to design propulsion strategies in the context of artificial micro-swimmers [45,86].

Applications
We now proceed to illustrate the capabilities of our solver with three different applications. We consider first a static problem in which self-contact, bending and twist give rise to the classic out-of-plane configurations denoted as plectonemes [40], while the addition of stretching and shearing produces a different type of experimentally observed solutions, known as solenoids [40]. Then we turn our attention to two dynamic biophysical problems in which an active filament interacts with a solid and a liquid environment, exhibiting qualitatively different optimal biolocomotion strategies. We emphasize here that one of the main points of these biolocomotion studies, besides verifying the versatility and prediction capabilities of our solver in biophysical settings, is to demonstrate its practical use in an inverse (optimization) design process, which typically represents a demanding stress test for any simulation algorithm, as unusual parameter sets or conditions are explored.

Plectonemes and solenoids
When an inextensible rod is clamped at one end and twisted a sufficiently large number of times at the other end, it becomes unstable, coils up and generates a characteristic structure known as a plectoneme [87]. While this behaviour has been well characterized both theoretically and experimentally [87], its analogue for highly extensible filaments has been largely ignored [40]. In particular, for large extensional and twisting strains qualitatively different solutions arise, such as those corresponding to tightly packed solenoidal structures [40] whose properties are as yet poorly understood, and whose importance has only recently been recognized in the context of soft robotics and artificial muscles [41][42][43][44].
Given the broad scope of our computational framework, we can now study the formation of both solenoids and plectonemes. As illustrated in figure 9a, a soft rod of Young's modulus E = 10 6 Pa is clamped at one end, and subject to an axial load F, while also being twisted R times at the other end. As experimentally and theoretically observed for F = 0, i.e. in the absence of stretching (L/L ≈ 1), plectonemes are generated (figure 9b). When the load F is increased so that the elongation of the rod approaches L/L ≈ 1.15, solenoids arise as predicted in [40] and illustrated in figure 9c. This test case, therefore, shows the ability of our solver to qualitatively capture different instability mechanisms, driven by the competition between the different modes of deformation of the rod. We leave the details of the explanation of the phase diagram for the formation of plectonemes, solenoids and intermediate structures [40] for a later study.

Slithering
The mechanics of slithering locomotion typical of snakes has been extensively investigated experimentally [83,88,89], theoretically [82,90,91] and computationally [92,93]. While biological experiments have provided quantitative insights, theoretical and computational models have been instrumental to characterize qualitatively the working principles underlying snake locomotion. Although these models implement different levels of realism, they generally rely on a number of key simplifications. Typically, theoretical models assume planar deformations [82] and/or disregard mechanics by prescribing body kinematics [91]. Computational models offer a more realistic representation, but they have mostly been developed for and tailored to robotic applications [92,93]. For example, snakes are often modelled as a relatively small set of hinges and/or springs representing pointwise localized actuators that connect contiguous rigid segments. Therefore, they do not account  for the continuum nature of elastic body mechanics and biological muscular activity. Moreover, in robot replicas the critical feature of friction anisotropy is commonly achieved through the use of wheels [89]. As a consequence computational models often assume only two sources of anisotropy, in the tangential and lateral direction with respect to the body. This is in contrast with biological experiments [83] that highlight the importance of all three sources of anisotropy, namely forward, backward and lateral.
Our approach complements those previous attempts by accounting for physical and biological effects within a continuum framework (equations (2.12)-(2.15)). In this section, we demonstrate the qualitative and quantitative capabilities of our solver by reverse engineering optimal slithering gaits that maximize forward speed.
We consider a soft filament of unit length actuated via a planar travelling torque wave of muscular activity in the direction perpendicular to the ground. The interaction with the substrate is characterized by the ratios μ f k , as experimentally observed for juvenile Pueblan milk snakes on a moderately rough surface [83]. The value of the friction coefficient μ f k is set so that the ratio between inertial and friction forces captured by the Froude number is Fr = (L/T 2 m )/(μ f k g) = 0.1, as measured for these snakes [83]. In the spirit of [53,56,57], we wish to identify the fastest gaits by optimizing the filament muscular activity. The torque wave generated by the snake is parametrized according to §4.1.1 and is characterized by N m = 6 control points and a unit oscillation period T m , so that overall we optimize for five parameters, four of which are responsible for the torque profile along the rod (β 1 , β 2 , β 3 , β 4 ), while the last one represents the wavenumber 2π/λ m (see §4.1.1).
These parameters are left free to evolve from an initial zero value, guided by an automated optimization procedure that identifies the optimal values that maximize the snake's forward average speed v fwd max over one activation cycle T m . The algorithm of choice is the Covariance Matrix Adaptation-Evolution Strategy [94,95] (CMA-ES) which has proved effective in a range of biophysical and engineering problems, from the optimization of swimming gaits [53], morphologies [56,57] and collective dynamics [55] to the identification of aircraft alleviation schemes [96] or virus traffic mechanisms [52]. The CMA-ES is a stochastic optimization algorithm that samples generations of p parameter vectors from a multivariate Gaussian distribution N . Here each parameter vector represents a muscular activation instance, and every generation entails the evaluation of p = 60 different gaits. The covariance matrix of the distribution N is then adapted based on successful past gaits, chosen according to their corresponding cost function value f = v fwd max , until convergence to the optimum. The course of the optimization is reported in figure 10  v fwd max 0.6, consistent with experimental evidence [97]. Moreover, CMA-ES finds that the optimal wavelength is λ m L (figure 10g), again consistent with biological observations [83,98]. Thus, this value of wavelength strikes a balance between thrust production and drag minimization within the mechanical constraints of the system. We note that a rigorous characterization of slithering locomotion would require the knowledge of a number of biologically relevant parameters (Young's and shear moduli of muscular tissue, maximum attainable torques, etc.) and environmental conditions (terrain asperities, presence of pegs, etc.) and goes beyond the scope of the present work. Nevertheless, this study illustrates the robustness, quantitative accuracy and suitability of our solver for the characterization of biolocomotion phenomena.

Swimming
We finally turn to apply the inverse design approach outlined in the previous section to the problem of swimming at low Reynolds numbers where viscous forces dominate inertial effects. We maintain the exact same set-up as in the slithering case, while we change the environment from a solid substrate to a viscous fluid. The flow-filament interaction is then modelled via slender-body theory, as illustrated in §4.1.5.
Once again we inverse design planar optimal gaits for forward average speed f = v fwd max within one activation cycle T m , by employing the same muscular activity parametrization as for slithering. To verify a posteriori the biological relevance of the identified optimal solution, we consider the case of the sea urchin spermatozoon Echinus esculentus [99], which swims by means of helical or planar waves travelling along its flagella of length L s 40 µm. The gait corresponding to planar swimming is characterized by kinematic undulations of wavelength λ s < L s and frequency f s 2.8. At Re 10 −4 the spermatozoon attains the scaled velocity v s = U s /(f s L s ) 0.08 ± 0.03, where U s is the dimensional cruise speed [99]. Although this gait may not be the absolute optimal planar locomotion pattern, the fact that it is replicated in a large number of organisms [85] suggests that it captures some effective features that we expect to qualitatively recover via our numerical optimization. The course of the optimization is reported in figure 11 together with the kinematic details of the identified fastest gait. As can be noted in figure 11e,f, the forward average scaled speed and wavelength approach v fwd max 0.055 and λ m 0.38L are qualitatively and quantitatively consistent with experimental evidence [99].
As in the previous section, we note that a rigorous characterization of optimal swimming at low Reynolds numbers would require the knowledge of a number of biologically relevant parameters and environmental conditions, a careful comparison with a rich body of literature [33,35,36,85], and goes beyond the scope of the present work. Nevertheless, this and the previous study illustrate how the interplay between filament mechanics and the surrounding environment crucially affects propulsive gaits, as is biologically evident and automatically recovered via our numerical inverse design approach.

Conclusions
We have presented a versatile implementation of the Cosserat rod model for the simulation of soft filaments deforming in three-dimensional space. These filaments at any given cross section can experience all deformation modes, namely normal and orthonormal bending, twisting, stretching and shearing. Particular emphasis is placed in realistically accounting for substantial axial extensions and shear strains, targeting emerging soft robotic applications. The solver also handles self-contact, muscular activity, solid boundaries, isotropic and anisotropic friction as well as hydrodynamics. The outcome is a simple algorithm suitable to tackle a plethora of engineering and biophysical phenomena.
We verified this against a battery of benchmark problems entailing different physical aspects and boundary conditions, and used it to solve a range of forward and inverse problems spanning mechanics and biophysics. In particular, we studied the formation of solenoids and plectonemes and the evolutionary optimization of terrestrial limbless locomotion and swimming.
Our results demonstrate the utility of a modelling approach based on Cosserat rods in a wide range of settings that involve active and passive soft filaments, broadening the scope of previous studies. We particularly emphasize that our robust computational approach allows us to approach the solution of an inverse problem by coupling our method with evolutionary strategies successfully, without the problems associated with the (often unforeseen) variety of candidate solutions produced throughout the process. Ongoing work involves its coupling to realistic high Reynolds number flow solvers [100], its integration with sensory feedback models for the characterization of locomotory neural circuitry [80], and its use for the rational design of biohybrid robots [45,101].
Data accessibility. Our software implementing the study cases presented in the manuscript and in the appendix is publicly available at https://github.com/mattialab/elastica.git. Acknowledgements. We thank Nicholas Charles for a careful reading of and comments on the paper.

Appendix A. Lagrangian governing equations
The following shows the conversion of governing equations for unstretchable filaments from the laboratory to material frame of reference. Note that in the Governing equations and Numerical method sections of the main text, all terms are scaled appropriately by the local dilatation to account for stretching.
The time and space derivatives of the centreline r(s, t) are associated with the velocity v and tangent field t v = ∂r ∂t and t = ∂r ∂s , ( A 1 ) with |t| = 1, since s is the current arc-length. Similarly, time and space derivatives of the material frame Q, due to the orthonormality of the directors, are associated by definition with the angular velocity ω and generalized curvature κ vectors, so that ∂d j and where the equivalences ∂ t Q T · Q = ω × (·) and ∂ s Q T · Q = κ × (·) hold. These kinematic equations combined with the linear and angular momentum balance laws at a cross section [102] yield the governing equations for the Cosserat rod and where ρ is the constant material density, A is the cross-sectional area in its current state (so that ρAv is the linear momentum line density), h(s, t) is the angular momentum line density, n(s, t) and τ (s, t) are the internal force and torque resultants, and f and c are external body force and torque line densities. The internal force n and torque τ resultants depend on the geometric and material properties of the filament. This dependence is embedded via the material or constitutive laws, which must be frame invariant, rendering their definition in a Lagrangian setting most natural. Moreover, we note that the angular momentum line density can be readily expressed in the material frame as Qh = h L = ρIω L , where the tensor I is the second area moment of inertia which, assuming circular cross sections, takes We can use these definitions to establish relationships between time and space derivatives of the laboratory and material frames. and Because the second moment of inertia is defined in the material frame, the governing equation for the angular momentum is also most naturally expressed in the material frame (I in the laboratory frame is a function of space and time rendering its use cumbersome, while I L is a constant and, in our isotropic case, diagonal matrix in the material frame). We can use equations (A 16) and (A 17) to convert all of the terms in equation (A 7) to the material frame. Multiplying by Q and noting that Iω = Q T I L ω L yields our final governing equation for the angular momentum, expressed in the material frame Aside from the balance of angular momentum, converting equations (A 5) and (A 6) to the material frame simply requires a direct application of the definition of the material frame to the RHS quantities.

Appendix B. Axial stretch comparison with Neo-Hookean and Mooney-Rivlin models
We consider the simple example of a beam subject to axial stretching and compare the load response predicted by our approach (that combines linear stress-strain constitutive laws and cross-sectional area rescaling) with Neo-Hookean [60] and Mooney-Rivlin models [61,62] for hyperelastic materials. We assume an initial lengthL = 100 mm and cross sectionÂ = 10 mm 2 . As the beam is stretched the cross-sectional area decreases. In one-dimensional, the Neo-Hookean model reads: where α is the stress, e = L/L is the elongation factor (related to the strain σ = e − 1), C 10 is a material property to be determined experimentally and that in the limit of small strains (e → 1) is related to the Young's modulus E, as above indicated. In one-dimensional, the Mooney-Rivlin model reads: lim e→1 E = 6(C 10 + C 01 ), where C 10 and C 01 are material properties to be determined experimentally and that in the limit of small strains (e → 1) are related to the Young's modulus E as indicated above. In our approach the load reads Assuming that the beam is made of rubber with E = 10 7 Pa, and that in the Mooney-Rivlin model C 10 = C 01 , we obtain the load responses illustrated in figure 12. As can be noted our approach reasonably approximates more advanced hyperelastic models within a 30% deformation range.

Appendix C. Time discretization
To advance the discrete system of equations (3.6)-(3.9), we choose the energy conserving, second-order position Verlet time integration scheme. This allows us to write a full iteration from time t to t + δt as and The time integrator is structured in three blocks: a first half step position update (equations (C 1) and (C 2)), followed by the evaluation of local accelerations (equations (C 4) and (C 5)), and finally, the second half step position update (equations (C 5) and (C 6)). Therefore, the second-order position Verlet scheme entails only one right-hand side evaluation of equations (3.8) and (3.9), the most computationally expensive operation. We also emphasize that for the numerical integration of equations (C 2) and (C 6) the Rodrigues formula is employed, implying the direct use of ω L instead of ω = Q T ω L . By foregoing an implicit integration scheme we can quickly incorporate a number of additional physical effects and soft constraints. This algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity.

Appendix D. Discrete operators
The discrete operators h : We note that h and A h operate on a set of N vectors and return N + 1 vectors of differences or convert integrated quantities to point-wise, respectively, consistent with equations (3.8) and (3.9).

Appendix E. Derivation of the vertical displacement of a Timoshenko cantilever beam
We briefly derive here the analytical expression of the vertical displacement y of equation (3.14) for the cantilever beam problem of figure 5a. To do so we make use of the free body diagram of figure 13, and of the constitutive relations of table 1. Moreover, we disregard axial extension and assume small deformations so that the coordinate x (along the direction k) is approximated by the arc-length s.
Recalling that the bending strain is the space derivative of the bending angle ψ (figure 13), we may write where M is the bending moment, E the Young's modulus, and I 1 is the second area moment of inertia about the axis j = k × i (figure 5a). The shear angle θ, as illustrated in figure 13, is the difference between the bending angle and the slope of the centreline, so that where V is the shear force, A is the cross-sectional area, G is the shear modulus and α c is the Timoshenko factor [58].

. Euler buckling instability
Euler buckling involves a single straight isotropic, inextensible and unshearable rod subject to an axial load F, as depicted in figure 14a. The critical axial load F c that the rod can withstand before bending can be expressed analytically [3] as where E is the modulus of elasticity of the rod, I = I 1 = I 2 is the area moment of inertia, L is the length and b is a constant which depends on the boundary conditions. If both ends are fixed in space but free to rotate then b = 1, thus where α = EI is the bending stiffness. We test our solver against the above analytical solution by simulating the time evolution of an inextensible and unshearable rod. In figure 14a, we show the results of our computations in the limit when both ends are free to rotate and all their spatial degrees of freedom are fixed except for one which allows the top end to slide vertically. The inextensible and unshearable conditions are enforced numerically by setting the entries of the matrix S to be much larger than those of B (details in figure 14). We explore the phase spaces F-α and F-L and determine the probability of a randomly perturbed rod to undergo a buckling event, characterized by the bending energy (defined in the main text as a quadratic functional) exceeding the small threshold value E B > E th . As can be seen in figure 14b,c, the obtained probability landscapes exhibit a sharp transition in agreement with the exact solution of equation (F 2).

F.2. Mitchell buckling instability
The Euler buckling benchmark allows us to test the capability of our solver to capture the onset of an instability relative to the bending modes. Next, we consider the Mitchell buckling instability that simultaneously accounts for both bend and twist. When the ends of an isotropic, inextensible and unshearable filament are joined together, the resulting configuration is a planar circular shape. If the two ends are twisted by a given angle and glued together, a circular shape with uniform twist density is obtained. When the total twist Φ, i.e. the integrated twist line density along the filament, exceeds a critical value Φ c , the rod buckles out of plane. This critical total twist can be expressed analytically [103] as where α and β are, respectively, the bending and twist rigidities. We explore the phase space F-(β/α) and determine the probability of a randomly perturbed rod to undergo a buckling event, characterized by the translational energy (defined in the main text as a quadratic functional) exceeding a small threshold value E T > E th . As can be seen in figure 15, our results show a sharp transition in agreement with the exact solution of equation (F 3).

F.3. Twist forced vibrations in a stretched rod
Next, we examine the problem of twisting waves generated in a rod axially stretched, as illustrated in figure 16a. The coupling between stretching and the other dynamic modes tests the rescaling in terms of the dilatation factor e of equations (2.14) and (2.15).   The rod is clamped at one end, stretched by the quantity L = (e − 1)L, and forced to angularly vibrate about the axial direction by applying at the free end the couple A v sin(2π f v t), where A v and f v are the corresponding forcing amplitude and frequency. In the case of no internal dissipation and constant circular sections, the induced standing angular wave θ(s, t) can be determined analytically. The dynamics of twisting yields the wave equation for the angular displacement [104] where c s = √ G/ρ is the shear wave velocity, G is the shear modulus and ρ is the density. By assuming a solution of the form θ (s, t) = φ(s) sin(2π f v t), and by applying the boundary conditions φ(0) = 0 and d s φ(0) = C v /(Î 3 G) with C v twisting torque andÎ 3 second area moment of inertia about the axial direction, we can solve equation (F 4) obtaining As can be seen in figure 16b, our numerical method recovers the derived analytical solution for the twist angular displacement along the stretched rod. Moreover, the solver consistently exhibits a second-order convergence in time and space given the error metric = θ − θ n , where θ n represents our numerical solutions in the limit of refinement (figure 16c).

Appendix G. Further friction validations G.1. Validation of rolling friction model
Here, we validate the friction model introduced in §4.1.5 on three test cases that can be analytically characterized. In all benchmarks, we consider a rigid, unshearable, inextensible straight rod of mass m, length L, radius r, and axial mass second moment of inertia J. The rod interacts with a surface characterized by static and kinetic friction coefficients μ s and μ k , thus assuming isotropic friction.
In the first test, the rod is initially at rest on a plane inclined of the angle α s , as depicted in figure 17a. Owing to the gravitational acceleration g the rod starts rolling or slipping, depending on α s , down the plane. The linear v and angular ω velocities of the filament, and therefore the corresponding translational E T = mv 2 /2 and rotational E R = Jω 2 /2 energies, can be analytically determined.
By recalling equation (4.12), the force necessary to ensure rolling without slip takes the form F noslip = −F /3 = −mg sin(α s )/3. Given the maximum static friction force F s = μ s F ⊥ = μ s mg cos(α s ) in the case |F noslip | ≤ |F s |, we have a = (F + F noslip )/m. Therefore, at the time T after releasing the rod, linear and angular velocities read, respectively, v = aT and ω =ωT = (a/r)T, expressing the no slip kinematic constraint between linear a and angularω accelerations. On the contrary, if |F noslip | > |F s | the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time T we have a = (F − μ k F ⊥ )/m, v = aT, ω =ωT = (μ k F ⊥ r/J)T. Therefore, translational and rotational energies as a function of α s finally read As can be noted in figure 17b, our numerical approach faithfully reproduces the derived analytical solution, accurately capturing the discontinuity at the transition between pure rolling and slipping. The second test case of figure 17c,d entails a rod set in motion by the external couple C s on a horizontal plane. Depending on the magnitude of the load the filament exhibits pure rolling or slipping motion. By recalling equation (4.12), the force necessary to ensure rolling without slip takes the form F noslip = 2C s /(3r). Given the maximum static friction force F s = μ s F ⊥ = μ s mg, in the case of |F noslip | ≤ |F s | we have a = F noslip /m. Therefore, at the time T after releasing the rod, linear and angular velocities read, respectively, v = aT and ω =ωT = (a/r)T, expressing the no slip kinematic constraint between linear a and angularω accelerations. On the contrary, if |F noslip | > |F s | the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time T we have a = μ k F ⊥ /m, v = aT, ω =ωT = J −1 (C s − μ k F ⊥ r)T. Therefore, translational and rotational energies as a function of C s finally read As can be noted in figure 17d, again our numerical approach faithfully reproduces the derived analytical solution, accurately capturing the discontinuity at the transition between pure rolling and slipping.
In the last test case, we consider a rod of initial velocity V s in the direction parallel to a horizontal plane and perpendicular to the filament axis (figure 17e), and we vary the relative mass moment of inertia ratio λ = 2J/(mr 2 ). Initially, the rod slips on the surface and gradually slows down, due to the effect of the kinetic friction, to a point for which linear a and angular velocity ω meet the kinematic constraint v eq = rω eq characteristic of pure rolling motion. The frictional force F induces the torque M = rF, so that over a period of time T we have v = V s − FT/m and ω = rFT/J. By enforcing the no slip kinematic constraint V s − FT/m = r 2 FT/J, we have that FT = V s /(r 2 /J + 1/m). This allows us to directly compute v eq and ω eq   Figure 18. Anisotropic static and kinetic longitudinal friction. A rod initially at rest on a horizontal plane is pulled longitudinally with force F. The system is characterized by anisotropic friction so that resistance forces in the forward direction F f long are larger than in the backward direction F b long . The plot illustrates the behaviour of the total rod's translational energy E T as a function of the force F applied. Blue and solid black lines correspond to, respectively, simulated and analytical E T for a rod pulled forward. Orange and dashed black lines correspond to, respectively, simulated and analytical E T for a rod pulled backward. Settings: length L = 1 m, radius r = 0.025 m, mass m = 1 kg, Young's modulus E = 10 5 Pa, shear modulus G = 2E/3 Pa, shear/stretch matrix S = 10 4 × 1 1 1 N, bend/twist matrix B = diag(EI 1 , EI 2 , GI 3 ) Nm 2 , dissipation constant γ = 10 −6 kg (ms) −1 , gravity g = 9.81 m s −2 , static and kinetic forward friction coefficients ( G 1 ) As depicted in figure 17f, our numerical approach accurately reproduces the predicted energies as a function of λ.

G.2. Validation of anisotropic longitudinal friction model
After validating our rolling friction model, we turn to its longitudinal counterpart. Here, we consider a straight, rigid, inextensible and unshearable rod which is axially pulled or pushed with force F for a fixed period of time T ( figure 18). Anisotropy is captured through the forward μ f s , μ f k and backward μ b s , μ b k static and kinetic coefficients. For F · t ≥ 0, the rod translational energy takes the form (|F| − μ f k mg) 2 if |F| > |μ f s mg| , while for F · t < 0 we have As can be seen in figure 18, our numerical method reproduces the above theoretical predictions.

G.3. Isotropic versus anisotropic friction driven locomotion
In this section, we illustrate the effect of symmetry breaking via anisotropic friction in a system constituted by a soft filament interacting via surface friction with a solid substrate. If we consider isotropic friction and a specular muscular activation pattern by setting the control values β 1 = β 4 and β 2 = β 3 and the wave number 2π/λ m = 0 (see §4.1.2), then the system is symmetric (up to inertial effects that can be neglected for small Froude numbers). Therefore, over any activation cycle the snake centre of mass does not move. On the contrary, if we introduce anisotropy the snake will be able to slowly move (the capability to effectively move is impaired by the absence of the travelling gait component since 2π/λ m = 0).
As can be observed in figure 19, this prediction is captured by our numerical scheme, which accurately resolves the physical mechanisms at play. Moreover, this test shows once again how our methodology is robust in terms of numerical noise as no spurious displacements or rotations are generated in the symmetric case ( figure 19a,b).

Appendix H. Representative timings
The proposed algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity: by foregoing an implicit integration scheme we can incorporate a number of additional physical effects and soft constraints, even though this may come at the expense of computational efficiency. Indeed, for large Young's or shear moduli or for very thin rods, the system of equations (2.4)-(2.7) in the main text might become stiff, so that small timesteps must be employed to ensure stability. Although we have not derived any rigorous CFL condition, throughout this work we employed the empirical relation dt = a d , with a ∼ 10 −2 s m −1 , and found it reliable in preventing numerical instabilities ( figure 20).

Representative timings
Optimal snake Optimal swimmer  Figure 20. Timings for various representative simulations. The time-to-solution values reported correspond to an average over 10 measurements for the same study case.
Despite the potential stiffness of this model, since the computational cost per timestep scales linearly with the number of discretization elements n (the model is one-dimensional), we could carry out all our computation on a single MacBook Pro equipped with a 2.3 GHz dual-core Intel Core i5, Turbo Boost up to 3.6 GHz, with 64 MB of eDRAM. Because of the condition dt = a d , a smaller number of elements implies larger timesteps (d = L/n). Therefore, the time-to-solution may scale between linearly and quadratically, depending on how the timestep is set. This can be noted by inspecting Euler and helical buckling studies reported in figure 20. In the Euler buckling simulations, the timestep is set to be the same for both case 1 and case 2, and indeed by halving n, the time-to-solution and the time per timestep reduce by a factor of 2. In the helical buckling study, the timestep was adaptively selected according to dt = a d , so that, when n is halved, the time-to-solution reduces by a factor of 4, while the time per timestep reduces by a factor of 2.