A study of snake-like locomotion through the analysis of a flexible robot model

We examine the problem of snake-like locomotion by studying a system consisting of a planar inextensible elastic rod with adjustable spontaneous curvature, which provides an internal actuation mechanism that mimics muscular action in a snake. Using a Cosserat model, we derive the equations of motion in two special cases: one in which the rod can only move along a prescribed curve, and one in which the rod is constrained to slide longitudinally without slipping laterally, but the path is not fixed a priori (free-path case). The second setting is inspired by undulatory locomotion of snakes on flat surfaces. The presence of constraints leads in both cases to non-standard boundary conditions that allow us to close and solve the equations of motion. The kinematics and dynamics of the system can be recovered from a one-dimensional equation, without any restrictive assumption on the followed trajectory or the actuation. We derive explicit formulae highlighting the role of spontaneous curvature in providing the driving force (and the steering, in the free-path case) needed for locomotion. We also provide analytical solutions for a special class of serpentine motions, which enable us to discuss the connection between observed trajectories, internal actuation and forces exchanged with the environment.

We examine the problem of snake-like locomotion by studying a system consisting of a planar inextensible elastic rod with adjustable spontaneous curvature, which provides an internal actuation mechanism that mimics muscular action in a snake. Using a Cosserat model, we derive the equations of motion in two special cases: one in which the rod can only move along a prescribed curve, and one in which the rod is constrained to slide longitudinally without slipping laterally, but the path is not fixed a priori (free-path case). The second setting is inspired by undulatory locomotion of snakes on flat surfaces. The presence of constraints leads in both cases to non-standard boundary conditions that allow us to close and solve the equations of motion. The kinematics and dynamics of the system can be recovered from a one-dimensional equation, without any restrictive assumption on the followed trajectory or the actuation. We derive explicit formulae highlighting the role of spontaneous curvature in providing the driving force (and the steering, in the free-path case) needed for locomotion. We also provide analytical solutions for a special class of serpentine motions, which enable us to discuss the connection between observed trajectories, internal actuation and forces exchanged with the environment.  robotics. This is a new and recent paradigm in robotic science [1,2], whereby inspiration is sought from nature to endow robots with new capabilities in terms of dexterity (e.g. the manipulation abilities of an elephant trunk or of an octopus arm) and adaptability (e.g. the ability of snakes to handle unexpected interactions with unstructured environments and move successfully on uneven terrains by adapting their gait to ground properties that change from place to place in an unpredictable way). The way snakes move has been the subject of seminal works by Gray [3,4]; see also [5,6]. In these early studies, Gray described the mechanics underlying snake locomotion inside closely fitting channels and on a surface in the presence of external push-points. Subsequently, muscular activity as well as forces transmitted by snakes to arrays of pegs among which they move have been measured [7,8]. Further early theoretical studies can be found in the Russian literature (e.g. [9][10][11][12], and the references quoted therein). More recently, focus has turned to the importance of frictional anisotropy between snakes' ventral skin and flat surfaces on which they move, stimulating both experimental and theoretical research [13][14][15][16][17]. In fact, it is well established that equality of friction coefficients in longitudinal and lateral directions leads to no net forward motion in undulatory locomotion (e.g. [15,16,18,19] for similar results in the closely related problem of undulatory swimming locomotion).

Introduction
The idea that frictional anisotropy plays a role in snake locomotion was put forward long ago in the engineering literature [5] and, most notably, by Hirose in his seminal work on robotic snakelike locomotion [20]. Hirose was among the first to realize the potential of biological inspiration in designing robots by studying snake-like locomotors and manipulators [20]. Technological advances in this field have led to the development of models for snake robots crafted with more and more jointed active segments, eventually leading to the use of continuum theories [21]. In some more recent contributions [22][23][24][25][26], Cosserat models are used for the mechanics of slender flexible robots, described as deformable rods.
Inspired by the literature on snake-like locomotion recalled above, in this paper we study a model system similar to the one used in [27] in the context of undulatory swimming, and consisting of a planar inextensible elastic rod that is able to control its spontaneous curvature. This is the curvature the rod would exhibit in the absence of external forces, which can be non-zero in the presence of internal actuation (see the sketch in figure 1b). Local control of this quantity provides an internal actuation mechanism that can be used to mimic muscular activity in biological undulatory locomotion. Indeed, by varying its spontaneous curvature α, the rod generates a distributed internal bending moment M a . The two quantities satisfy the simple relation M a = −EJα, where EJ is the bending stiffness of the rod (see (2.4)). Travelling waves of spontaneous curvature can put the system in motion when the environment exerts constraints or forces that prevent the rod being deformed everywhere according to its spontaneous curvature.
To show how control of spontaneous curvature in the presence of external constraints leads to locomotion, we use a Cosserat model and derive the equations of motion for two special cases: one in which the rod can only move along a prescribed curve (prescribed-path case), and one in which the rod is constrained to slide longitudinally without slipping laterally, but the path is not fixed a priori (free-path case). The first case corresponds to a rod confined in a channel with frictionless walls. The second case is inspired by the slithering motion of snakes that interact through anisotropic frictional forces with a flat surface on which they are free to move. Frictional resistance is typically larger in the lateral direction than in the longitudinal one. Our setting corresponds to the limiting case of infinite ratio between lateral and longitudinal friction coefficients, in which longitudinal sliding is allowed while lateral slipping is forbidden.
Our work is closely related to the approach presented in [17], which we extend in at least one major way. In fact, in [17] locomotion of an active rod with no lateral slipping along a free path is considered. The authors impose, however, periodicity of the solution (effectively considering a rod of infinite length) which leads to an incomplete system of equations. The system is then closed by postulating laws (closure relations, justified by experimental observations) on the lateral forces exerted on the ground surface. The novelty of our approach consists in solving the equations of motion in the case of a system of finite length, with no a priori assumptions either on the followed path, which can be non-periodic, or on the reactive forces imposing no lateral slipping. These both emerge as part of the solution of the problem, once a history of spontaneous curvatures is assigned. Closure of the equations is obtained by carefully considering edge-effects, which lead to non-standard boundary conditions. We derive in this way explicit formulae that enable us to explore in full generality the connection between observed motion, internal actuation and lateral forces exchanged with the environment. Moreover, we are able to solve inverse locomotion problems; namely, given a motion of the system that we want to observe, find an internal actuation that produces it.
Our main results are the following. We formulate direct and inverse locomotion problems (direct: find the motion produced by a given actuation history; inverse: find the actuation history required to produce a given motion), and show existence and uniqueness of the solution of direct problems, and non-uniqueness for the inverse ones. In the prescribed-path case, we reduce the dynamics of the system to a single ordinary differential equation for the tail-end coordinate (the only degree of freedom for an inextensible rod forced to slide along a given curve). This equation reveals clearly the mechanism by which a flexible rod can actively propel itself inside a channel, whenever the channel exhibits a variation of curvature along its track, and provides a quantitative framework to revisit some of the classical findings on snake motility by Gray. In the free-path case, we are again able to close the equation of motion and reduce the dynamics of the system to a single equation, this time an integro-differential equation for the tail-end coordinate. A particularly interesting outcome of our analysis is the emergence of an asymmetry in the mechanical boundary conditions at the (leading) head and the (trailing) tail. This is not only a mathematical subtlety, but it is also deeply grounded in the physics of the problem. While the tail follows the path traced by the preceding interior points, the head is free to veer laterally, 'creating' the path as the motion progresses. We show that the curvature of this newly created path is set by the time history of spontaneous curvatures at the leading head. Recognizing this steering role of the spontaneous curvature leads to a procedure to generate solutions for the freepath case from those of the prescribed-path case, based on modifying them near the leading head, in order to account for steering. Again, we provide explicit formulae to calculate the lateral forces transmitted to the ground surface.
The rest of the paper is organized as follows. In §2, we present our mathematical model of a flexible robot as an active rod, and formulate direct and inverse locomotion problems. In §3, we derive the governing equations and the appropriate boundary conditions for motion inside a channel with frictionless walls (prescribed-path), solve them in some simple geometries and discuss the physical implications of our results. In §4, we derive the governing equations and corresponding boundary conditions for the motion of an active rod sliding longitudinally without slipping laterally on a flat surface (free-path) and propose a class of analytical serpentine solutions. Possible connections of our results with observations made in the context of biological snake locomotion are briefly summarized in §5, while the existence and uniqueness of the solution of the equations of motion for the free-path case is proved in appendix A.

The flexible robot model
We consider a model consisting of a (long) chain of cross-shaped elements (figure 1b) linked together by ideal joints connected by deformable springs. We assume that each spring is able to actively change its rest length (the length at which the tension in the spring is zero). Following [22][23][24][25][26], we model this system through a continuous description based on the planar Cosserat rod theory.
A configuration of a Cosserat rod of reference length L on the plane is defined by a pair of vector-valued functions where b is a unit vector. The curve r describes the midline of the rod, while b characterizes the orientation of its deformed cross sections (figure 1a). As in [28], we introduce also the unit vector a := −e 3 × b, where e 3 is the unit vector normal to the plane. We then define the strain variables ν and η through the following decomposition along the moving orthonormal frame {a, b}: where the subscript s is used to denote the partial derivative with respect to the space variable. The function ν = ν(s, t) describes the stretch, while η = η(s, t) defines the shear strain. Finally, the bending strain μ := θ s is obtained through the scalar valued function θ(s, t) defined by where {e 1 , e 2 } is a fixed basis in the plane containing the rod. We consider our system as being made of an infinite number of elements such as those in figure 1b, each of them being of infinitesimal length, and assembled along the central curve r of the rod. As we assume them to be rigid, we impose the constraints that the rod is inextensible and unshearable: The ability of the robot to modify the equilibrium length of each of the connecting springs can be naturally modelled macroscopically by considering an elastic rod which can actively vary its spontaneous curvature; namely, the curvature the rod would exhibit in the absence of external loads. This is similar to what is done [27] in the context of swimming motility. We model this by introducing the elastic potential density where EJ is the bending stiffness of the rod. Note that if (2.2) hold, then r s always coincides with the unit vector a, and the bending strain μ(s, t) is equal to the curvature of the rod at the point r(s, t). Therefore, the function α in (2.3) can be viewed as a varying spontaneous curvature, which we assume to be freely controllable in order to set the robot in motion. The bending moment and can be seen as the sum of a passive elastic term EJθ s and of an active one M a := −EJα which can be varied at will by suitably tuning α. An active moment originating from muscular contraction is used in the model of snake locomotion in [17]. Along with the elastic potential we define the kinetic energy density where the subscript t denotes the partial derivative with respect to time, ρA is the linear mass density and ρJ is the linear moment of inertia. Finally, the Lagrangian density L of the system reads In the following sections, we will consider two types of locomotion problems arising from the interaction of prescribed spontaneous curvature and external constraints. The direct one can be formulated as follows: given a time history of spontaneous curvatures α(s, t), together with initial and boundary conditions, find the motion r(s, t) of the rod and the forces it exchanges with the environment. In the inverse one, the motion is prescribed, and we want to find a history α(s, t) that produces it, together with the corresponding forces. We will consider two types of external constraints and see that, in both cases, the direct problem has a unique solution while, for the inverse one, the solution is not unique. For studies of swimming locomotion problems conducted in a similar spirit, we refer the reader to [19,27,[29][30][31].
3. The case of prescribed path: sliding inside a channel The first problem we consider is motion along a prescribed path. We place our robot model inside a curved channel fitting exactly its body, and we assume that there are no friction forces exerted by the walls of the channel. We model such a setting by imposing the external (holonomic) constraint where the equation φ Γ = 0 defines (we assume, globally) the curve Γ , which we interpret as the central line of the channel. There is no loss of generality in assuming |∇φ Γ | = 1.

(a) Derivation of the equations of motion
We derive the equations of motion through Hamilton's principle, adding to (2.5) an external for all variations δr and δθ defined on [0, L] × [t 1 , t 2 ] and vanishing at its boundary. If we define n := Na + Hb, the Euler-Lagrange equations we obtain from (3.2) are where the bending moment M is defined in (2.4). These are the classical dynamical equations for a planar Cosserat rod (e.g. [28]) with an external force given, in our case, by the transversal reaction imposing constraint (3.1). We can suppose that our active rod is in frictional contact with the ground. The presence of a longitudinal frictional force per unit length is handled by simply adding F , where Sgn denotes the sign function, to the left-hand side of the first equation.
To close the equations of motion we use the principle of mechanical boundary conditions (PoMBC) [32]. We define generalized edge loads acting on the system by considering the rate at which work is expended at the edges in virtual motions compatible with the constraints, and assume that all generalized edge loads acting on the system are explicitly prescribed.
In view of (3.1), we have that where s 0 and s L are the curvilinear coordinates relative to Γ of the two ends of the rod, which we call generalized edge coordinates, and Θ is the angle between the tangent vector to Γ and e 1 , so that Using (3.4) to derive the expressions for r t and θ t at s = 0, L, we obtain where we used a 'dot' to denote the time derivative of the generalized coordinates and k is the curvature of Γ . The coefficients multiplying the generalized velocitiesṡ 0 (t) andṡ L (t) are the generalized edge loads which, by the PoMBC, have to be prescribed. As we suppose that no external edge forces are doing work on the system at either of the two ends, we enforce the condition P edge = 0 by setting such loads equal to zero. Finally, conditions (2.2) and (3.1) must be added to the equations of the system. As the active rod is assumed to be inextensible and unshearable, and its backbone curve r is forced inside the graph of Γ , the constrained system can be described with only one degree of freedom, namely the curvilinear coordinate relative to Γ of the first end of the robot model. Thus, and substituting these expressions in the equations of motion we obtain, accounting also for longitudinal friction, where k = k(s 0 (t) + s). As for the boundary conditions, they now read Summarizing, in order to solve the (direct) locomotion problem stated at the end of §2, we need to find the unknown functions N(s, t), H(s, t), f (s, t) and s 0 (t). The equations we have for this purpose are the three equations of motion (3.8)-(3.10), and the two boundary conditions (3.11). We see that, by integrating (3.8), a first-order ordinary differential equation (ODE) in the space variable s, we can derive one additional ODE (in the time variable) containing only the unknown s 0 (t), which completely determines the motion of the system. This ODE is given below as equation (3.13), or (3.14) in a simplified version. Once s 0 is known, we can use (3.10), (3.8) and (3.9), together with the boundary condition (3.11) holding at s = 0, to determine H, N and f , respectively.
We show now how to obtain the ODE for s 0 (t). If we substitute in (3.8), the expression of H given by (3.10), then integrating on the space variable, we have where m = L 0 ρA ds is the total mass of the rod, (a) Two configurations of the elastic rod inside a spiral channel: initial (light grey) and final (dark grey). A positive work W is necessary to vary the position of the end point from Γ (ξ 2 ) to Γ (ξ 1 ) and force the rod inside the channel. (b) Upon release, the first endpoint slides back from Γ (ξ 1 ) to Γ (ξ 2 ) and the rod exits the channel with velocity V.
If we now apply (3.11), we obtain the equation which, complemented with initial position and velocity, defines s 0 uniquely. The shear force H is now uniquely defined by (3.10), while Using all the expressions above, we can recover f from (3.9).
Let us now suppose that our active rod is stiff and slender enough, so that EJ, ρA ρJ. We can then neglect the terms containing ρJ in (3.13), obtaining the simplified equation (3.14) Equation (3.14) shows that the dynamics of the rod is reduced to that of a point particle of mass m subjected to a force given by the sum of three terms. The first one is a 'potential' force depending exclusively on the geometry of Γ , the second one is a friction term, while the third is an 'active' force which depends on the spontaneous curvature α.
The following examples illustrate the role played by these terms in the dynamics of the system.

(b) Spiral channel
Let us consider only the first term in the right-hand side of (3.14) by setting α, γ = 0. The system described in this case is a passive elastic rod with straight rest configuration (α = 0) placed inside a curved channel with frictionless walls and no frictional interaction with the ground (γ = 0). Observe that the only non-vanishing term in the right-hand side of (3.14) is the first one, which states that the driving force on the rod depends only on the curvature of the channel at the two ends of the body (this can be interpreted as a result of inextensibility). Moreover, the sign of this force is such that the rod is always pushed towards the region of smaller curvature. As an example, consider the case of a spiral-shaped channel where k(s) = K/s, with K > 0 (figure 2). Then (3.14) with α = 0 reads In order to thread the rod inside the spiral by varying the coordinate of the endpoint from ξ 2 to ξ 1 , we have to do positive work as we have to increase the curvature at every point of the body. If we then release the rod it will accelerate towards the exit and return back to ξ 2 with a positive velocity The system moves towards a 'straighter' configuration, decreasing its elastic energy and therefore increasing its kinetic energy. Similar problems of passive elastic rods sliding inside frictionless sleeves have been studied, both analytically and experimentally, in [33]. Let us suppose that α, γ = 0. The elastic rod can now vary its spontaneous curvature and it has to overcome a longitudinal frictional force to slide inside the spiral channel. The active force term in (3.14) can assume any value if we suppose that we have no restrictions in the choice of α. Thus, in particular, an active elastic rod can slide inside the spiral without need of external pushing. More generally, the system can achieve motion in a predetermined direction when placed inside any channel which does not present circular or straight sections of length greater than L. This last result is reminiscent of the theoretical and experimental findings of J. Gray in his study [3] of snake undulatory locomotion. Using an energy balance argument, he concludes that it is possible for a snake to slide inside a channel closely fitting its body only provided such a channel exhibits a variation of curvature along its track. He then shows experimentally that snakes are able to move in sinusoidal closely fitting channels, but motion in straight ones only occurs through a different gait (concertina), which is impossible if the width of the channel and of the snake body are comparable.
We consider now two concrete examples of an active rod propelling itself inside the spiral channel. We do this by solving an inverse locomotion problem: we prescribe the motion of the rod and then deduce two histories of spontaneous curvatures that produce it. In this way, we also show non-uniqueness of the inverse problem.
Suppose we want to find an activation that propels the rod inside the spiral at constant velocitẏ s 0 (t) = −V < 0. Equation (3.14) then reads where s 0 (t) = s in − Vt for some initial value s in for s 0 at t = 0. Set then an easy calculation shows that (3.  Estimating lateral forces associated with internally actuated conformational changes, and minimizing them, may be of interest in the field of minimally invasive interventional medicine, e.g. for concentric-tube continuum robots, also called active cannulas [26]. For small values of the geometric parameter ζ , the channel is close to a straight tube while, as ζ grows, it becomes wavier and wavier. The wavelength λ dictates how many turns the channel has per unit length. We want to find a history of spontaneous curvatures α(s, t) that produces motion along the sinusoidal channel (3.18) with constant longitudinal velocitẏ Assuming that the trailing edge of the active rod lies at the origin at t = 0, we must have s 0 (t) = Vt. If we also assume that L = 2π nλ, where n is a positive integer, then the potential term in equation (3.14) vanishes, and constant forward motion is realized only if the active force exactly matches the frictional one: We give two different solutions for the history of spontaneous curvatures α satisfying (3.20) (i.e. generating the same prescribed motion) by solving two constrained minimization problems. Among all α's such that (3.20) holds, we find the ones that minimize I bend (the bending energy) and I act (the activity), where To solve, for example, the second problem we consider the extended functional where q is the Lagrange multiplier enforcing (3.20). The spontaneous curvature α bend minimizing the bending energyÎ bend must then solve δÎ bend [α bend ; q] = 0, where the variation of the extended functional is taken with respect to α. A straightforward calculation gives where the second equality is obtained by plugging the expression for α bend in (3.20). More explicitly, using (3.19), we get From the equations of motion and the boundary conditions, taking again ρJ = 0, we then obtain Note that none of the external and internal forces depend on the bending stiffness EJ. This allows us to consider the rigid limit EJ → ∞, for which the observable motion and forces do not vary, whereas, on the other hand, α bend (s, t) → k(s + Vt). This limit case could be relevant for the steering of wheeled robots in which curvature control is achieved through internal motors. Let us find α act that minimizes I act by repeating the procedure above. We obtain in this case that the optimal α is proportional to the derivative of the channel's curvature k s , whereby internal actuation is concentrated around inflection points of the trajectory. This is reminiscent of patterns of muscular activity observed in snake undulatory locomotion [7,8] for internal and external forces generated by α act , in terms of the corresponding quantities we found for α bend . We observe that the two force fields differ by terms proportional to EJ, while they become indistinguishable when EJ → 0. We give here a graphical representation of the two solutions, using material parameters taken from the zoological literature. Based on [8], we set L = 1.3 m and γ = μ mg L −1 , where μ = 0.2 is the longitudinal friction coefficient, m = 0.8 kg and g is the gravitational acceleration constant. Following [15], we neglect the inertial terms in all the expressions, setting ρAV 2 = 0. As for the bending stiffness, we explore a range going from EJ = 10 −4 Nm 2 [34] to EJ = 10 −3 Nm 2 [17].
Results are shown in figure 4 for n = 3 and ζ = 18.5 m −1 , while we choose the largest value of EJ in order to emphasize the difference between the two solutions in terms of forces exerted on the channel walls. The force field f bend consistently displays maxima in magnitude near the inflection points of curvature. On the other hand, f act varies substantially during motion: at some times it displays local maxima at the points of maximal concavity and convexity, while at some other times maxima are located at the inflection points. Note that, at points of maximal concavity and convexity of the rod, the lateral force f is perpendicular to the horizontal axis (the average direction of motion) and does not contribute to propulsion. We will comment further on these features in the next sections.

The free-path case
We now turn to the case in which the path is not a priori known and study an active rod free to move on a flat surface through longitudinal sliding without lateral slipping. Accordingly, we impose the (non-holonomic) constraint where r ⊥ s = e 3 × r s . We denote by −f r ⊥ s the transversal reactive force per unit length (exerted by the ground on the rod) enforcing the no-slip condition, where f is the Lagrange multiplier associated with constraint (4.1). At the same time, we suppose that a frictional force F given by (3.3) acts in the longitudinal direction.
This choice for F relies on the simplifying assumption that frictional forces are uniform along the rod's body. Moreover, real systems such as snakes [13][14][15] or snake-like robots [20] cannot rely on transversal frictional reactions of arbitrary magnitude to prevent lateral slipping. Solutions of interest for a more realistic description of these systems can be considered; for instance, those for which the reactive force f imposing constraint (4.1) does not exceed a maximum value, which can be determined experimentally.

(a) Derivation of the equations of motion
We deduce the equations of motion through the Lagrange-d'Alembert principle, similar to what is done in [35,36]. The principle states that, in the presence of the dissipative force F , a solution (r, θ) that satisfies constraint (4.1) must solve We complement these equations with boundary conditions by relying again on the PoMBC. Let us consider a typical configuration of the robot in motion while subjected to the external constraint (4.1). We suppose that such a movement is directed head-first, where we denote the head as r(L, t) and the tail as r(0, t). As shown in figure 5a, an asymmetry between head and tail emerges. Because of (4.1), the tail position and director can change only by assuming the values previously taken at an adjacent internal point. We can therefore impose on s = 0 the same conditions we had in the channel case, namely where v 0 is the (only) generalized velocity at s = 0 and k(0, t) is the curvature of r evaluated in s = 0 at time t. As for the head, as the path is no longer predetermined, we have an extra degree of freedom. Condition (4.1) requires r t and r s to be collinear; therefore, this extra degree of freedom must come from the rotation of the director. We then impose where v L and ω L are the generalized velocities for the system at s = L. The work rate of the external edge forces is Thus, there are two generalized edge loads at s = L, namely the axial tension n · r s and the bending moment M, and one at s = 0, with the same expression it had in the channel case. We set the generalized loads equal to zero because, just like in the previous section, we suppose that no external edge force is doing work on the system. Alongside with the boundary conditions coming from the vanishing of the generalized edge loads, the system must be complemented with equations (2.2) and (4.1).
The non-holonomic constraint (4.1) compels the active rod to move within a curve in the plane, much like it was for the channel case in the previous section. This time, however, the path is not a priori determined but is created during the motion, and it is an unknown of the problem. In fact, constraints (2.2) and (4.1) lead to the existence of some function s 0 and some curve Γ , which have both to be determined, such that (3.7) holds. As the boundary conditions we derived hold only for head-first motions, we only consider solutions satisfyinġ The equations of motion written in components are formally identical to the ones derived for the channel case and EJ(k s − α s ) + H = ρJ(ks 0 (t) + k sṡ0 (t) 2 ), (4.7) but k = k(s 0 (t) + s) is no longer predetermined. On the other hand, the equations are closed through the boundary conditions obtained by setting the three generalized edge loads equal to zero N(0, t) + EJ(k(s 0 (t)) − α(0, t))k(s 0 (t)) = 0 and N(L, t) = 0, EJ(k(s 0 (t) + L) − α(L, t)) = 0, (4.8) together with the initial curvature k(s 0 (t 0 ) + s) = K(s), with s ∈ [0, L], and the initial values for s 0 andṡ 0 at t = t 0 . Such values must satisfy the compatibility relationṡ s 0 (t 0 ) > 0, K(L) = α(L, t 0 ) and K s (L)ṡ 0 (t 0 ) =α(L, t 0 ). (4.9) In order to solve the locomotion problem, we need to find the unknown functions N(s, t), H(s, t), f (s, t), together with s 0 (t) and k(s 0 (t) + s). The three equations of motion and the three boundary conditions (4.8) are sufficient to solve this problem uniquely. This leads to a unique solution also for r and θ , once the initial position r(0, t 0 ) and orientation θ(0, t 0 ) of the first end are prescribed, by integrating the equations θ s = k and r s = Γ as done, for example, in [15,16]. The detailed proof is provided in appendix A, and we only sketch here the heuristic argument behind it.
A key role is played by the third boundary condition in (4.8), coming from the vanishing of the bending moment at the leading edge. This latter condition, namely k(s 0 (t) + L) = α(L, t), (4.10) assigns a crucial role to the spontaneous curvature at the leading edge in determining the path followed by the system. Thus, the value of α at s = L operates as a 'steering wheel', while the internal values of the spontaneous curvature supply the active force for propulsion, as for the channel case. Let us see how s 0 and k can be determined. There is no loss of generality if we take t 0 = 0 and s 0 (0) = 0. On the other hand, let us assumeṡ 0 (t) > 0 for t ∈ [0, t * ) so that s 0 is invertible in the whole interval, and let us also assume that t * is small enough so that s 0 (t) < L for every t. Clearly, k(s) = K(s) is known for s ∈ [0, L] from the initial condition. For s > L, we can recover k from the history of spontaneous curvatures at the leading edge because each point of the path Γ (ξ ) with ξ > L generated between t 0 and t * is the location of the leading edge at some time s −1 0 (ξ − L) (figure 5b). Thus, setting we can recover k(s 0 (t) + s) from the initial conditions, the given α and the knowledge of s 0 . In turn, s 0 can be determined by substituting the expression for H given by (4.7) into (4.5) and integrating with respect to s. Using (4.8), we deduce (m + ρJQ(s 0 (t)))s 0 (t) = EJ 2 (k 2 (s 0 (t)) − k 2 (s 0 (t) where R and Q are given by (3.12). Moreover, using (4.11) and the change of variable s = ξ − s 0 (t), the last integral in (4.12) can be written as the sum The second summand in the last expression can be rewritten further, using the change of variable ξ = s 0 (τ ) + L, as where we have used the identity k s (s 0 (t) + L)ṡ 0 (t) =α(L, t) following from (4.10). Finally, observing that, in view of our assumption s 0 (t) < L, we have k(s 0 (t)) = K(s 0 (t)) and k(s 0 (t) + L) = α(L, t), it follows that R =ṡ 0 (t) 2 2 (α 2 (L, t) − K 2 (s 0 (t))) and Q = L s 0 (t)

Equation (4.12) is in fact
an integro-differential equation in s 0 alone which can be uniquely solved in terms of the data of the problem, as proved in appendix A. Just like in the channel case, once s 0 and k are known, the unknown functions H, N and f can be readily deduced from (4.7), (4.5) and (4.6), respectively.

(b) Serpentine solutions
In this section, we provide a class of explicit serpentine solutions for the free-path locomotion problem, by exploiting solutions constructed for the channel case. We obtain these solutions by solving an inverse locomotion problem, prescribing the motion first and then looking for a history of spontaneous curvatures α(s, t) that produces it.
Let us consider the sinusoidal path Γ given by (3.18) and assume that our active rod slides at constant longitudinal velocity V, so that s 0 (t) = Vt. As we did before, we set ρJ = 0 for simplicity. Following the arguments of §3c, we conclude that α must again solve (3.20). In addition, we must now also require the boundary condition (4.10), which assigns the steering role to α, to be satisfied. Note that none of the spontaneous curvatures we obtained in the channel case fulfils (4.10). However, as we show in the following, we can locally modify any α solving (3.20) so that (4.10) is also satisfied.
We focus below on the history of spontaneous curvatures α act given by (3.22), because it is the one that more closely resembles the typical muscular activity patterns observed in undulating snakes. If we consider a function α(s, t) = α act (s, t) +α(s, t) (4.13) with a 'steering term'α such that L 0α (s, t) cos s + Vt λ ds = 0, (4.14) then α satisfies both (3.20) and (4.10). With α having these properties, s 0 (t) = Vt becomes a solution for the equations of motion, and the expression for N, H and f can be deduced following the procedure of §3a. The termα in (4.13) satisfying (4.14) can be taken of the form where δ is an arbitrary constant, which can be set as small as we want, and p i (t) with i = 3, . . . , Q are coefficients explicitly depending on t and implicitly depending also on δ and all the dynamical parameters. These coefficients can be uniquely determined by imposing (4.14) and any other Q − 5 linearly independent relations between them (e.g. in the numerical experiment we are about to propose, we imposedα ss (L, t) = 0, which led to a smooth generated force field f concentrated near the head). Note that the functionα so defined is twice continuously differentiable. If we take δ small enough, then α differs from α act only in a small neighbourhood of the leading edge where the steering termα is non-zero. The reactive shear force and tension are now given by H(s, t) = H act (s, t) + EJα s (s, t) and while the force exerted on the ground reads, in this case, From the last equalities, it follows that, if δ is small, forces (external and internal) have the same values of the corresponding ones obtained in the channel case with the exception of a small region near the leading edge. We set δ/L = 0.25 and we give here two graphical comparisons of the same solution fitted with different parameters (figures 6 and 7).
In figures 6a and 7a, we take the same values we considered in §3c for all the parameters. When compared with that of figure 4b, this solution clearly shows the asymmetry that the steering termα generates in the activation and force pattens in the proximity of the head (leading edge). In figure 6b, we take the smaller value for the bending stiffness, EJ = 10 −4 Nm 2 . Note that this solution displays a similar force pattern to that of figure 4a (as expected from the formulae we derived in §3c), which is generated by a different spontaneous curvature history. Finally, in figure 7b, we consider an active rod with the same bending stiffness but moving with a less tortuous gait (smaller ζ ). Observe that, also in this case, we obtain an almost stationary force pattern, which is qualitatively similar to that of figure 6b.
Summarizing, we see a picture consistent with that of snake undulatory locomotion hypothesized in [17] (muscular activity and lateral forces both concentrated near the inflection points of the trajectory, where the propulsive effect of the lateral forces is largest because their component along the average direction of motion is largest) that emerges either automatically, for specific choices of material parameters (figure 6b) or through adjustment of the gait (figure 7b). Lateral forces near points of maximal and minimal convexity may also be ruled out by eliminating ground contact (by lifting portions of the body near those points), as is done in [15,16] and sometimes observed in undulating snakes.

Discussion
We have studied the motion of an active rod (a planar inextensible elastic rod of finite length with adjustable spontaneous curvature), arising from the interaction between external constraints and internal actuation by spontaneous curvature. Using Cosserat theory, we have formulated and solved both direct and inverse locomotion problems for two cases: one in which the system is forced to move along a prescribed path, and the other in which the path is not fixed a priori and the system slides along its tangential direction while subjected to lateral forces preventing lateral slipping. We have obtained a procedure to generate free-path solutions from solutions with prescribed-path, by recognizing the dual role (pushing and steering) played by spontaneous curvature in powering undulatory locomotion of the rod. Finally, we have obtained explicit analytic solutions and formulae that can be used to study the connections between observed motion, internal actuation and forces transmitted to the environment, and to explore how these connections are affected by the mechanical properties of the system (its bending stiffness).
Although our results hold for a (very specific) model system, it may be interesting to compare some of them with observations made in the context of undulatory locomotion of snakes. For this exercise to make sense, we are formulating the implicit assumption that our mechanism of internal actuation by spontaneous curvature can provide a reasonable proxy for muscular actuation, and that the free-path motion of the organism we are considering does not cause lateral slipping, but only involves longitudinal sliding (as is sometimes observed).
The first example is equation (3.14), which provides a compact summary of some classical observations on snake locomotion by Gray [3,4]. Undulatory locomotion in closely fitting channels is possible only if the channel presents a variation of curvature along its track. The formula explains the mechanism by which spontaneous curvature can provide the driving force for locomotion inside a tightly fitting channel, and our analysis delivers formulae to calculate the lateral forces exerted on the channel walls. It would be interesting to compare these with experimental measurements.
A second example is the observation that, among various possible actuation strategies producing the same prescribed motion, the one minimizing actuation effort (as measured by the integral norm of spontaneous curvature) is proportional to the arc-length derivative of the curvature of the trajectory. This means that local actuation is maximal at the inflection points of the trajectory, and zero at points of maximal and minimal curvature. This is closely reminiscent of the typical pattern of muscular actuation emerging from experimental measurements on snakes [7,8], and it would be interesting to explore further the reasons behind this analogy.
Finally, our analysis suggests that the connection between observed motions, internal actuation and transmitted forces may be strongly affected by the passive mechanical properties of the system, such as its bending stiffness. The conceptual picture of snake undulatory locomotion in which both muscular activity and lateral forces are concentrated near the inflection points of the trajectory, previously theorized in [17], can emerge either automatically, for specific choices of material parameters, or through the adjustment of the gait Understanding the mechanisms that control gait selection and, in particular, whether there are optimality criteria explaining it in biological organisms, and whether some of them may be useful for the engineering of artificial devices, represent interesting challenges for future work (see, however, [16,17], for results in this direction). Adding some important ingredients, currently not present in our model, may prove necessary. One example is some form of active local control of the frictional interactions between body and ground, as is done in [15,16]. Moreover, when considering real snakes' behaviour it is natural to speculate that muscular activity may be, at least to some extent, a reaction to external stimuli (the forces exerted by the ground on the snake), thereby creating an interplay between the two dynamical variables. It would be interesting to study how our model could be extended to account for such feedback mechanisms. All these questions will require further study.
Data accessibility. All computations were performed in Mathematica, v. 10.0 (Wolfram Research, Inc.). This article does not report any primary data.
Authors' contributions. The authors have contributed equally to this work. Competing interests. The authors declare they have no competing interests. Funding. This work was supported by the European Research Council through the ERC advanced grant no.
( A 1 ) The general case in which s 0 (t) > L may also occur (i.e. the trailing edge is no longer contained in the image of the initial configuration, figure 5b) can be handled by applying the following simple, yet technical, argument. As we will show, a local solution s 0 for (A 1) exists and is unique once a positive initial velocityṡ 0 (t 0 ) is given, by requiring that s 0 (t) ≤ L. If a local solution s 0 of (A 1) exists and is unique then, for all given initial conditions, there exists only one solution with a maximal interval of definition. For such maximal solutions we can have either of two cases. In the first, the maximal interval of existence of s 0 withṡ 0 > 0 is of the type [0, t * ) and, for every t in the interval, s 0 (t) ≤ L holds. If that occurs, the only solution of the free-path problem is defined in the time interval [0, t * ), the curvature k can be derived through (4.11) while all the other unknowns can be deduced by the same procedure we employed in the channel case in §3a. In the second case, a solution s 0 of (A 1) satisfyingṡ 0 > 0 and s 0 (t) ≤ L can only be defined in a maximal domain of the type [0, t * ]. For a solution of this kind we must have s 0 (t * ) = L, by maximality. In this last case we can still define k through (4.11) as we did before. Then we can reapply all the arguments of §4a finding an equation of the type (A 1) for a new variable s * 0 with new initial conditions for the free-path locomotion problem, namely s * 0 (t * ) = s 0 (t * ),ṡ * 0 (t * ) =ṡ 0 (t * ) and the new initial curvature K * (s) = k(L + s) with s ∈ [0, L]. After that we are able to solve the new integro-differential problem uniquely for s * 0 , recover the value for all the unknowns, and then glue the solutions together. We repeat this procedure until we reach eventually a maximal domain of existence for the general solution.
The existence and uniqueness of free-path locomotion solutions then follows from the local existence and uniqueness of solutions of (A 1) with the extra requirementsṡ 0 (t) > 0 and s 0 (t) ≤ L. This can be proved using standard contraction mapping arguments. The result holds under the For small enough T the operator C is a contraction from B M,T x 0 into itself, therefore it has only one fixed point. This proves local existence and uniqueness for the extended version of (A 2). If we take x 0 = (0, y 0 ) with y 0 > 0 then, restricting the domain of existence to an interval [0, T * ) if necessary, we have by continuity x(t) ≤ L andẋ(t) = y(t) > 0 for every t ∈ [0, T * ), hence obtaining the unique solution to the original (non-extended) problem.