## Abstract

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.

## 1. Introduction

Snake locomotion has fascinated natural scientists for a long time. More recently, it has become a topic of great interest as one of the key examples of soft bioinspired 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–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–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).

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 snake-like 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–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 1*b*). 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 free-path 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.

## 2. The flexible robot model

We consider a model consisting of a (long) chain of cross-shaped elements (figure 1*b*) 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–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
**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 1*a*). 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**}:
*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
**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 1*b*, 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:
**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 resulting from (2.3) is
*θ*_{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
*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 *N*=*N*(*s*,*t*) and *H*=*H*(*s*,*t*) are the reactive internal forces (axial tension and shear force, respectively) enforcing constraints (2.2).

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–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
*ϕ*_{Γ}=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 reactive potential −*fϕ*_{Γ}(**r**), where *f*=*f*(*s*,*t*) is the Lagrange multiplier enforcing (3.1). A solution (**r**,*θ*) must satisfy
*δ***r** and *δθ* defined on [0,*L*]×[*t*_{1},*t*_{2}] and vanishing at its boundary. If we define **n**:=*N***a**+*H***b**, the Euler–Lagrange equations we obtain from (3.2) are
*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
**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
*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

*P*

_{edge}of the edge loads as

**r**

_{t}and

*θ*

_{t}at

*s*=0,

*L*, we obtain

*k*is the curvature of

**. The coefficients multiplying the**

*Γ**generalized velocities*

*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,**

*Γ**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
*s*_{0} uniquely. The shear force *H* is now uniquely defined by (3.10), while
*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
*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
*ξ*_{2} to *ξ*_{1}, we have to do positive work
*ξ*_{2} with a positive velocity

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 velocity *s*_{0}(*t*)=*s*_{in}−*V* *t* for some initial value *s*_{in} for *s*_{0} at *t*=0. Set
*α*=*α*_{1} and *α*=*α*_{2}. Moreover, if the history of spontaneous curvatures *α* is given by either *α*_{1} or *α*_{2}, then *s*_{0}(*t*)=*s*_{in}−*V* *t* becomes automatically a solution for the equations of motion (3.8)–(3.10), and *N*, *H* and *f* can be explicitly written following the procedure illustrated in the previous section. We denote by *f*_{1} and *f*_{2} the lateral forces exerted by the active rod under the actuations *α*_{1} and *α*_{2}, respectively. In figure 3, the two solutions are illustrated with *L*=0.5 m, *K*=8 m^{−1}, *γ*^{′′}=0.3, EJ=10^{−3} Nm^{2} and where we set *ρAV* ^{2}*k*_{s}=0, thereby ignoring inertial effects. An interesting consequence of this last assumption is that *f*_{1}=*γ*^{′′}(*K*+2/*K*) is constant, as figure 3*a* shows. By contrast, *f*_{2} is not constant and it even changes sign (figure 3*b*).

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].

### (c) Sinusoidal channel

We address in this section an inverse locomotion problem for a sinusoidal channel (figure 4) meandering around the horizontal axis
*ζ*, 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 velocity
*t*=0, we must have *s*_{0}(*t*)=*V* *t*. 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:
*α* 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
*q* is the Lagrange multiplier enforcing (3.20). The spontaneous curvature *α*_{bend} minimizing the bending energy *α*. A straightforward calculation gives
*α*_{bend} in (3.20). More explicitly, using (3.19), we get
*ρJ*=0, we then obtain
*α*_{bend}(*s*,*t*)→*k*(*s*+*V* *t*). 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]. More in detail,
*q* given again by (3.21). In order to compare the two solutions, we write
*α*_{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.

## 4. 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
*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–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
*δ***r** and *δθ* that vanish at the boundary of [0,*L*]×[*t*_{1},*t*_{2}], while *δ***r** also satisfy
*δθ* must vanish, while the coefficient relative to *δ***r** must take the form *f*=*f*(*s*,*t*) is the unknown Lagrange multiplier enforcing constraint (4.1). The equations of motion then read

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 5*a*, 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
*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
*v*_{L} and *ω*_{L} are the generalized velocities for the system at *s*=*L*. The work rate of the external edge forces is
*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 satisfying

*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

*k*(

*s*

_{0}(

*t*

_{0})+

*s*)=

*K*(

*s*), with

*s*∈[0,

*L*], and the initial values for

*s*

_{0}and

*t*=

*t*

_{0}. Such values must satisfy the compatibility relations

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
*α* 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 *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

*b*). Thus, setting

*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

*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

*ξ*=

*s*

_{0}(

*τ*)+

*L*, as

*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

*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*)=*V* *t*. 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
*α* satisfies both (3.20) and (4.10). With *α* having these properties, *s*_{0}(*t*)=*V* *t* 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 *δ* 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 *f* concentrated near the head). Note that the function

If we take *δ* small enough, then *α* differs from *α*_{act} only in a small neighbourhood of the leading edge where the steering term *δ* 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 6*a* and 7*a*, we take the same values we considered in §3c for all the parameters. When compared with that of figure 4*b*, this solution clearly shows the asymmetry that the steering term *b*, 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 4*a* (as expected from the formulae we derived in §3c), which is generated by a different spontaneous curvature history. Finally, in figure 7*b*, 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 6*b*.

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 6*b*) or through adjustment of the gait (figure 7*b*). 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.

## 5. 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. 340685-MicroMotility.

## Acknowledgements

This study was started after an inspiring lecture given by Professor D. Bigoni at the ‘Material Theories’ workshop held in 2013 at the Mathematisches Forschungsinstitut Oberwolfach. We thank Professor F. L. Chernousko for pointing out to us important references in the Russian literature. A.D.S. thanks F. Renda for useful discussions.

## Appendix A. Existence and uniqueness of free-path solutions

For the sake of simplicity, we take *ρJ*=0. The following arguments can be easily adjusted for the general case.

Suppose we have a solution of the free-path locomotion problem (4.5)–(4.7) satisfying the boundary conditions (4.8) and the extra requirements (4.4) and (4.9), where *K*(*s*)=*k*(*s*_{0}(*t*_{0})+*s*) for *s*∈[0,*L*]. Again, there is no loss of generality taking *t*_{0}=0 and *s*_{0}(*t*_{0})=0. Let us first suppose, by restricting its domain of definition if necessary, that *s*_{0} is defined on a time interval [0,*t**) such that *s*_{0}(*t*)≤*L* for every *t*∈[0,*t**). The arguments used in §4a show that, within this time interval, the equation for *s*_{0} in terms of the data of the problem reads
*s*_{0}(*t*)>*L* may also occur (i.e. the trailing edge is no longer contained in the image of the initial configuration, figure 5*b*) 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 *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 *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 *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 *K**(*s*)=*k*(*L*+*s*) with *s*∈[0,*L*]. After that we are able to solve the new integro-differential problem uniquely for

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 *s*_{0}(*t*)≤*L*. This can be proved using standard contraction mapping arguments. The result holds under the very reasonable assumption of *α* and *K* being differentiable and uniformly bounded together with their derivatives.

Observe that we can recast (A 1) into a set of integro-differential equations of the form
**x** solves (A 2) if and only if *s*_{0}(*t*):=*x*(*t*) solves (A 1). We first extend *K* and *α* outside **x**(0)=**x**_{0} and no extra assumption on the solutions besides differentiability. The problem can be easily proved to be equivalent to that of the existence of a fixed-point for the operator *C* defined as
*C* to the space *t*↦**x**(*t*) defined on the interval *t*∈[0,*T*] and such that
*K* and *α* lead to the existence of two constants *M*_{G} and *M*_{H} such that
*L*_{G} and *L*_{H} such that
*τ* and *t* and for every *T* the operator *C* is a contraction from **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*∈[0,*T**), hence obtaining the unique solution to the original (non-extended) problem.

- Received January 26, 2015.
- Accepted November 4, 2015.

- © 2015 The Authors.