Front. Robot. AI Frontiers in Robotics and AI Front. Robot. AI 2296-9144 Frontiers Media S.A. 10.3389/frobt.2018.00099 Robotics and AI Original Research Peristaltic Waves as Optimal Gaits in Metameric Bio-Inspired Robots Agostinelli Daniele 1 Alouges François 2 DeSimone Antonio 1 3 * 1International School for Advanced Studies (SISSA), Trieste, Italy 2Centre de Mathématiques Appliquées, École Polytéchnique, Université Paris-Saclay, Paris, France 3The BioRobotics Institute, Sant'Anna School for Advanced Studies, Pisa, Italy

Edited by: Matteo Cianchetti, Scuola Sant'Anna di Studi Avanzati, Italy

Reviewed by: Shinya Aoi, Kyoto University, Japan; Barbara Mazzolai, Fondazione Istituto Italiano di Technologia, Italy

 *Correspondence: Antonio DeSimone desimone@sissa.it

This article was submitted to Soft Robotics, a section of the journal Frontiers in Robotics and AI

 05 09 2018 2018 5 99 15 04 2018 31 07 2018 Copyright © 2018 Agostinelli, Alouges and DeSimone. 2018 Agostinelli, Alouges and DeSimone

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Peristalsis, i.e., a motion pattern arising from the propagation of muscle contraction and expansion waves along the body, is a common locomotion strategy for limbless animals. Mimicking peristalsis in bio-inspired robots has attracted considerable attention in the literature. It has recently been observed that maximal velocity in a metameric earthworm-like robot is achieved by actuating the segments using a “phase coordination” principle. This paper shows that, in fact, peristalsis (which requires not only phase coordination, but also that all segments oscillate at same frequency and amplitude) emerges from optimization principles. More precisely, basing our analysis on the assumption of small deformations, we show that peristaltic waves provide the optimal actuation solution in the ideal case of a periodic infinite system, and that this is approximately true, modulo edge effects, for the real, finite length system. Therefore, this paper confirms the effectiveness of mimicking peristalsis in bio-inspired robots, at least in the small-deformation regime. Further research will be required to test the effectiveness of this strategy if large deformations are allowed.

 crawling motility lumbricus terrestris peristalsis self-propulsion metameric robots biomimetic robots soft robotics optimization MicroMotility - 340685 European Research Council10.13039/501100000781 1. Introduction

The study of self-propelled locomotors exploiting friction-induced traction as a result of body shape changes, is gaining attention because of the variety of physical systems which take advantage of such a locomotion strategy. One motivation is the desire to understand biological phenomena, such as cell migration on or within solid substrates, matrices, and tissues (Alberts et al., 2002). Another motivation is the attempt to replicate these mechanisms in robotics with the idea that biomimetic constructs may outperform traditional ones when confronted with unstructured and unpredictable environments.

In particular, robotic locomotion research has recently considered crawling and burrowing animals (e.g., earthworms, snakes, and caterpillars), whence an increasing number of research projects on bio-inspired metameric (soft) robots (Menciassi et al., 2006; Wang et al., 2009; Boxerbaum et al., 2012; Daltorio et al., 2013; Fang et al., 2015; Nemitz et al., 2016; Umedachi et al., 2016; Ge et al., 2017). As a matter of fact, many species such as earthworms, caterpillars, sea cucumbers and snails move using peristalsis which is a locomotion mechanism consisting of a series of wave-like muscle relaxation and contraction which propagate along the body (Quillin, 1999). One of the most studied biological species is Lumbricus terrestris (commonly known as nightcrawler) which is a kind of earthworm which uses peristalsis both for surface crawling and for burrowing. Each of its metameres (body segments) is endowed with longitudinal and circular muscles and can regulate frictional forces thanks to microscopic bristles called setae (Quillin, 1999). Understanding how relatively simple organisms are able to attain peristalsis and to which extent coordination is regulated by either the nervous system or spontaneous reflexes, are questions addressed by researchers for about a century and are still drawing attention (Garrey and Moore, 1915; Gray and Lissmann, 1938; Gardner, 1976; Quillin, 1999).

In the field of robotics, peristalsis has been mostly mimicked by a priori assignment of “gaits” defined by a few scalar parameters. Optimization of locomotion performances with respect to variations of these scalar parameters has been studied. Fang et al. (2015) consider harmonic deformations with a single, fixed, (time) frequency and amplitude, and determine the phase patterns of actuation maximizing the average velocity. Optimization leads to phase coordination, in the form of a pattern which is close to the identical-phase-difference (IPD) pattern corresponding to peristalsis. However, no rigorous proof of the connection between peristaltic waves and optimal actuation is given and, more importantly, the basic hypothesis of harmonic oscillations with a single fixed time frequency and amplitude is taken as an a priori assumption.

This paper aims to provide a deeper understanding of harmonic oscillations and peristalsis as result of an optimization problem rather than an a priori hypothesis. Indeed, we prove that - in the regime of small deformations - peristalsis is a symmetry property of the solution to an optimization problem. Symmetry of the solution comes from symmetry properties of operators in the equations governing the optimization problem, which are, in turn, the signature of geometric symmetries of the physical system.

The rest of the paper is organized in three main sections: material and methods, results and discussion. The first one, inspired by nightcrawlers' retractable setae, introduces a velocity-force law which is able to describe, for limit values of a single scalar parameter, the linear case (Newtonian) as well as the case of “free slip–perfect grip.” This model is presented in both continuous and discrete version and in the latter case we address some related optimal control problems. The second section illustrates the behavior of the continuous model by means of two examples and shows how peristalsis emerges, up to edge-effects, from optimization problems in the discrete framework. The last section presents a comparison with the work presented by Fang et al. (2015), states the main conclusions and provides directions for the future.

 2. Materials and methods 2.1. Continuous self-propelled 1D crawlers Model description and kinematics

Consider a 1D crawler moving along a straight line and assume its reference configuration is the segment

[X1:=0,X2:=L],

cf. Figure 1. We use the same formalism introduced by DeSimone and Tatone (2012) and DeSimone et al. (2013) so that X is the coordinate in the reference configuration while x(t) denotes the coordinate along the crawler's body in the current configuration (at time t). In particular x(t) is the image of X through the current transformation χ(·, t) which can be written in terms of the current distance s(·, t) from the left end, i.e.,

x(t)=χ(X,t):=x1(t)+s(X,t)∈[x1(t),x2(t)]

where x1(t): = χ(X1, t) and x2(t): = χ(X2, t). By definition, s(0, t)≡0 for all t and we assume that

s′(X,t):=∂s(X,t)∂X>0 ∀X,t

in order to guarantee the monotonicity of χ(·, t) at any time t.

Kinematics of a continuous 1D crawler: reference (A) and current (B) configurations.

In what follows a prime will denote the derivative with respect to X while a superscript dot will denote the derivative with respect to t.

We define the displacement

ua(X,t):=χ(X,t)−X

so that

χ′(X,t):=∂χ(X,t)∂X=1+ua′(X,t)=1+ϵ(X,t)

where we have implicitly defined the strain

ϵ(X,t):=∂ua∂X,

in terms of which condition (1) reads

ϵ(X,t)>−1 ∀X,t.

Finally, notice that, since the material (or Lagrangian) velocity is

χ˙(X,t)=x˙1(t)+s˙(X,t),

the spatial (or Eulerian) velocity is given by

v(x,t)=χ˙(X,t)|X=χ−1(x,t). Equations of motion

Throughout this section we deal with the motility problem, namely, given a history of strain ϵ(X, t), the aim is to find x1(t) which determines the dynamics of the one-dimensional crawler.

Friction laws

The force at the interface between substrate and crawler is modelled through a force-velocity relationship. In particular, we write the density per unit current length of the tangential component of the friction force at time t, f(x, t), as a function of the Eulerian velocity v(x, t).

Several models for the resistance forces are conceivable such as a

Newtonian model, i.e., a linear viscous law f(x,t):=−μv(x,t)

where μ>0 is a friction (or viscosity) coefficient;

or a more general “p-model” f(x,t):=−μgp(ϵ(X,t))v(x,t)     where gp(ϵ):=(11+ϵ)p

for p∈[0, +∞). Parameter p in our force law (4) allows us to investigate different types of frictional behaviors. For p = 0, we obtain a force per unit current length that is a linear function of velocity alone, which reduces to the Newtonian model (3). For p>0, we obtain a friction law that is sensitive to the state of elongation of the segment, with force per unit length higher or lower than that of the Newtonian case depending on whether the element is contracted (λ < 1 or ϵ < 0) or extended (λ>1 or ϵ>0). In the limit p → ∞, this produces an idealized model for friction in which no force opposes slip when the segment is extended (free slip), while the segment can withstand any tangential force without sliding (perfect grip) when it is contracted. We call this idealized model “free slip–perfect grip.” Figure 2 displays the graphs of gp(ϵ) around ϵ = 0 for different values of p. In fact, our model is a continuous analog of the discrete model proposed by Fang et al. (2015) to mimic the behavior of earthworms' setae, which protrude when the body is axially contracted, resulting in an increment of the resistance (Edwards et al., 2018).

Function gp(ϵ) governing the friction law (4) for selected values of parameter p.

In what follows, we will use the p-model (4).

 Force balance

The total friction is obtained by integrating the force per unit current length on the whole current domain, i.e.,

Ff(t)=∫x1x2f(x,t)dx=∫0Lfref(X,t)dX

where

fref(X,t):=f(x1(t)+s(X,t),t)s′(X,t).

Then, by neglecting inertia, the force balance yields

0=Ff(t)+Fe(t)=−μ∫0L(11+ϵ(X,t))p(x˙1(t)+s˙(X,t))s′(X,t)dX+Fe(t)=[−μ∫0L(11+ϵ(X,t))ps′(X,t)dX]x˙1    −μ∫0L(11+ϵ(X,t))ps˙(X,t)s′(X,t)dX+Fe(t)=[−μ∫0L(1+ϵ(X,t))1−pdX]x˙1    −μ∫0L(1+ϵ(X,t))1−ps˙(X,t)dX+Fe(t).

The square bracket multiplying ẋ1(t) in the formula above is the drag for rigid motion at unit speed and fixed shape ϵ(X, t), while Fe(t) is an external force which, for instance, can take into account the gravity force acting on a crawler on an inclined plane. Solving for ẋ1(t), we obtain

x˙1(t)=−∫0L(1+ϵ(X,t))1−ps˙(X,t)dX∫0L(1+ϵ(X,t))1−pdX              +Fe(t)μ∫0L(1+ϵ(X,t))1−pdX

which, in the case of zero external forces, is independent of μ.

Notice that, once the initial position x1(0) the strain ϵ(X, t) and the external force Fe(t) are provided, the whole dynamics x1(t) can be determined by integrating (5). Indeed, since

s(X,t)=s(0,t)+∫0Xs′(Y,t)dY                =∫0Xϵ(Y,t)dY+X

one has

s˙(X,t)=∫0Xϵ˙(Y,t)dY

and hence the right hand side of (5) is known once ϵ(X, t) is specified.

 2.2. Discrete self-propelled 1D crawlers

We now move to a discrete model, directly inspired by studies on annelid worms (Quillin, 1999) and metameric robots (Menciassi et al., 2006; Daltorio et al., 2013; Fang et al., 2015).

Model description and kinematics

We model the crawler's body as made up of N segments of same length L in the reference configuration (cf. Figure 3)

(Xn−1:=(n−1)L, Xn:=nL)

for n = 1, …, N. Let x0(t) denote the current position of the left edge; any point X in the reference domain is mapped to a point x(t) in the current domain through the map

x(t):=χ(X,t)=x0(t)+s(X,t)=X+ua(X,t)

whence the definition of strain

ϵ(X,t):=ua′(X,t)=s′(X,t)−1 .

Kinematics of a discrete 1D crawler consisting of N identical segments of reference length L. (A) Reference configuration. (B) Current configuration.

Furthermore, we assume that each segment can be contracted or expanded according to a constant stretch so that the overall strain results to be a piecewise constant function of X (at any fixed time t), i.e.,

ϵ(X,t):={ϵ1(t) if X∈[X0,X1) (segment 1)⋮ ϵN(t) if  X∈(XN−1,XN](segment N)

Consequently its time-derivative, ϵ°(·,t), is piecewise constant and hence, from (7), ṡ(·, t) is piecewise affine, i.e.,

s˙(X,t)=L∑j=1n−1ϵ˙j(t)+[X−(n−1)L]ϵ˙n(t)

for X∈[Xn−1, Xn]. Note that in this new framework the monotonicity condition (2) reads

ϵn(t)>−1 for all t and for n=1,…,N

which is the only constraint for an admissible history of strains, the datum of our motility problem.

 Equations of motion

Analogously to the previous case, the force balance yields

0=[−μ∫0NL(1+ϵ(X,t))1−pdX]x˙0(t)        −μ∫0NL(1+ϵ(X,t))1−ps˙(X,t)dX+Fe(t)

where, in view of (9),

∫0NL(1+ϵ(X,t))1−pdX=L∑j=1N(1+ϵj)1−p

and, in view of (10),

∫0NL(1+ϵ(X,t))1−ps˙(X,t)dX    =∑n=1N(1+ϵn)1−p·∫(n−1)LnL[​(X−(n−1)L)ϵ˙n(t)+L∑k=1n−1ϵ˙k(t)​]dX    =L22∑n=1N(1+ϵn)1−p[ϵ˙n(t)+2∑k=1n−1ϵ˙k(t)]    =L22∑j=1N[(1+ϵj)1−p+2∑m=j+1N(1+ϵm)1−p]ϵ˙j(t).

Solving for ẋ0(t) we obtain

x˙0(t)=[∑j=1N(1+ϵj)1−p]−1{​Fe(t)μL                −L2∑j=1N[​​(1+ϵj)1−p+2∑m=j+1N(1+ϵm)1−p​​]​​ϵ˙j(t)​}

which can be rewritten in the following vectorial form

x˙0(t)=∑n=1Nvn(ϵ(t))ϵ˙n(t)+q(ϵ(t))Fe(t)            =v(ϵ(t))·ϵ.(t)+q(ϵ(t))Fe(t)

where

vn(ϵ)=−L2(1+ϵn)1−p+2∑m=n+1N(1+ϵm)1−p∑j=1N(1+ϵj)1−p,q(ϵ(t))=[μL∑j=1N(1+ϵj(t))1−p]−1,v(ϵ)=[v1(ϵ)⋮vN(ϵ)] and ϵ(t)=[ϵ1(t)⋮ϵN(t)].

Equation (11) fully describes the dynamics once x0(0) and ϵ(X, t) are provided. In particular, the displacement after T time units is given by

Δx0=x0(T)−x0(0)         =∫0Tv(ϵ(t))·ϵ.(t)dt+∫0Tq(ϵ(t))Fe(t)dt. Relative displacement

We can rewrite everything in terms of a displacement relative to x0(t), defined as

u(X,t):=ua(X,t)−x0(t)                 =s(X,t)−X                 =∫0Xϵ(Y,t)dY

in order to describe the displacement in a coordinate system which is “co-moving” with the left end x0(t).

In the discrete framework, the relative displacement turns out to be a piecewise affine function of X (at any fixed time t): if X∈[Xn−1, Xn],

u(X,t)=L∑j=1n−1ϵj(t)+[X−(n−1)L]ϵn(t).

Setting

un(t):=u(Xn,t) for  n=1,…,N,

we have

ϵn=un−un−1L for  n=1,…,N,

where u0≡0. Equivalently

ϵ(t)=Ju(t)

where ϵ=(ϵ1,ϵ2,…,ϵN)T, u=(u1,u2,…,uN)T and

J:=1L[1   −11   ⋱⋱   −11] . Optimal control problems

In this section we address the problem of maximizing the net displacement Δx0 among periodic shape changes ϵ(X, t) with the same given energy cost.

We now describe the optimization problems with quadratic energy in the non-linear case first, and then in the small-deformation regime, for which general results can be established. We assume Fe≡0.

Feasible region

We assume that the shape function ϵ(t) is a C2 function defined from ℝ to ℝN.

In addition, we require ϵ(·) to be a time-periodic function. Finally we restrict our search to shape functions with a given cost per period, i.e., E[ϵ,ϵ°]=c, where the energy functional is assumed to be of the following quadratic form (in both ϵ and ϵ°)

E[ϵ,ϵ.]:=∫0TAϵ·ϵdt+∫0TBϵ.·ϵ.dt

where 𝔸 and 𝔹 are symmetric and positive definite N-dimensional matrices. Overall, the feasible region is

S:={ϵ∈C2(ℝ,ℝN) | ϵ(0)=ϵ(T)∧ |E[ϵ,ϵ.]=c}. Optimization problem

The general (non-linear) optimization problem is

maxϵ∈S Δx0[ϵ]:=∫0Tv(ϵ(t))·ϵ˙(t)dt

which is an isoperimetric problem (e.g., Van Brunt, 2004) involving N dependent variables ϵn. The corresponding Euler-Lagrange equations lead to a second order non-linear system of ODEs, i.e., for n = 1, …, N,

ddt∂ℱ∂ϵ˙n(t,ϵ,ϵ˙)−∂ℱ∂ϵn(t,ϵ,ϵ˙)=0

where F(t,ϵ,ϵ°):=v(ϵ)·ϵ°-λ(Aϵ·ϵ+Bϵ°·ϵ°) (here λ denotes a Lagrange multiplier).

 The small-deformation regime

We can focus on the small-deformation regime by expanding the objective function at the leading orders (about ϵ = 0), i.e.,

Δx0=∫0Tv(ϵ)·ϵ˙dt         =∫0T(v(0)+vϵ(0)ϵ+o(ϵ))·ϵ˙ dt         ≃v(0)·(ϵ(T)−ϵ(0))+∫0Tvϵ(0)ϵ·ϵ˙ dt         =∫0Tvϵ(0)ϵ·ϵ˙ dt

and, integrating by parts,

∫0Tvϵ(0)ϵ·ϵ˙ dt=[vϵ(0)ϵ·ϵ]0T−∫0Tvϵ(0)ϵ˙·ϵ dt                             =∫0T−ϵ˙·vϵT(0)ϵ dt

whence

Δx0≃12[∫0Tvϵ(0)ϵ·ϵ˙dt+∫0T−vϵT(0)ϵ·ϵ˙dt]          =∫0T(vϵ(0)−vϵT(0)2)ϵ·ϵ˙dt=:V[ϵ,ϵ˙].

In particular, it can be proved that the (skew-symmetric Toeplitz) matrix skw(vϵ(0)) =:V depends only on N, L and p. Indeed,

{V}ij={L(p−1)i−j+N2N2if i<j0if i=j−L(p−1)j−i+N2N2if i>j

(see Appendix A in Supplementary Material).

Therefore, in the regime of small deformations, problem (18) can be replaced by the following linear problem

maxϵ∈S V[ϵ,ϵ˙]:=∫0Tϵ˙·Vϵdt.

The corresponding Euler-Lagrange equations

ddt∂ℒ∂ϵ˙(t,ϵ,ϵ˙)−∂ℒ∂ϵ(t,ϵ,ϵ˙)=0

where

ℒ(t,ϵ,ϵ˙):=Vϵ·ϵ˙−λ(Aϵ·ϵ+Bϵ˙·ϵ˙),

lead to the following system of second order linear ODEs

Vϵ˙=λ(Bϵ¨−Aϵ) .

In general, a solution to (21) might be difficult to determine due to the complexity of finding a common diagonalization of 𝔸 and 𝔹. However, following the procedure adopted by Wiezel et al. (2018), we can solve this problem when one of the two operators is null, say 𝔸≡0 (resp. 𝔹≡0), and the other one, 𝔹 (resp. 𝔸), is symmetric, positive definite and such that the eigenspaces associated with the maximum-modulus eigenvalues of 𝔹-12𝕍𝔹-12 (resp. 𝔸-12𝕍𝔸-12) have dimension 1. Indeed, as shown in sections 1 and 2 in Appendix B (Supplementary Material), it turns out that

for A = 0 and B symmetric and positive definite, up to a constant, a solution of (20) must be of the form ϵ(t)=−Tπℜ(αie2πiTte)

where α∈ℂ\{0} is a constant such that ||α||=c2T and e=(e1,e2,…,eN)T∈ℂN\{0} is a suitable constant vector depending only on A and V.

for 𝔸 symmetric and positive definite and 𝔹 = 0, a solution of (20) with ϵ of unitary time frequency must be of the form ϵ(t)=2ℜ(αe2πiTte)

where α∈ℂ\{0} is a constant such that ||α||=c2T and e=(e1,e2,…,eN)T∈ℂN\{0} is a suitable constant vector depending only on 𝔹 and 𝕍.

Both (22) and (23) have the form

ϵ(t)=ℜ(α^e2πiTte)

i.e., they are circles in the plane (ℜ(e), ℑ(e)), regardless the number of links. Moreover, using the polar representations

α^=ϱaeiϑa and en=ϱneiϑn ∀n,

we get, for n = 1, …, N,

ϵn(t)=ϱaϱnℜ(ei(2πTt+ϑa+ϑn))          =ϱaϱnsin(2πTt+ϑa+ϑn+π2),

i.e., the optimal gait depends only on the 2N+2 parameters {ϑn}n, {ϱn}n, ϑa and ϱa. Admittedly, since α is a constant with fixed modulus and free argument, we can always assume ϑa+π2=0, i.e.,

ϵn(t)=ϱaϱnsin(2πTt+ϑn),

thus reducing the number of parameters to 2N+1.

The problem for 𝔸 = 0 and 𝔹 = 𝕀N is essentially equivalent to the one for 𝔸 = 𝕀N and 𝔹 = 0 (provided that unitary time frequency of ϵ is prescribed). Indeed, if (25) is a solution to

                                maxϵ∈S V[ϵ,ϵ˙]S:={ϵ∈C2 | ϵ(0)=ϵ(T)∧∫0T||ϵ||ℝN2dt=1}

then it is a solution also to

                                      maxϵ∈S V[ϵ,ϵ˙]S:={ϵ∈C2 | ϵ(0)=ϵ(T)∫0T∧  ∫0T||ϵ˙||ℝN2dt=(2πT)2}

and vice versa.

In general the two problems, 𝔸 = 0 with 𝔹 symmetric positive definite and 𝔹 = 0 with 𝔸 symmetric positive definite, are not equivalent. In fact, constraining the norm induced by one operator does not determine the norm induced by the other one, but only provides a bound. Indeed, denoting by λmin(·) and λmax(·) the minimum and maximum eigenvalue respectively, observe that, for ϵ(t) like (25),

∫0TAϵ·ϵdt≥λmin(A)∫0Tϵ·ϵdt                     =λmin(A)(T2π)2∫0Tϵ˙·ϵ˙dt                     ≥(T2π)2λmin(A)λmax(B)∫0TBϵ˙·ϵ˙dt

and, analogously,

∫0TBϵ˙·ϵ˙dt≥λmin(B)∫0Tϵ˙·ϵ˙dt                     =λmin(B)(2πT)2∫0Tϵ·ϵdt                      ≥(2πT)2λmin(B)λmax(A)∫0TAϵ·ϵdt . 3. Results 3.1. Continuous model: two examples of contraction waves

In this section we discuss two examples of contraction waves to illustrate the behavior of the p-model. We show that the parameter p determines the kind of motion: for p < 1 the motion is prograde (i.e., motion in the same direction as the one of the waves) while for p>1 the model reproduces an earthworm-like retrograde motion (i.e., motion in the opposite direction as the one of the waves).

For simplicity, in the following examples we neglect external forces, i.e., Fe≡0.

Smooth contraction wave

Consider a smooth traveling contraction wave by prescribing the strain along the body of the crawler as

ϵ(X,t):=ϵ0cos(2πL(X−ct))

or equivalently, in terms of the stretch,

s′(X,t)=1+ϵ(X,t)=1+ϵ0cos(2πL(X−ct))

where ϵ0 is the wave amplitude, L is the reference length of the crawler and c is a parameter which modulates time frequency and it is assumed to be strictly positive, i.e., the wave travels toward the right. By integrating over space,

s(X,t)=s(0,t)+∫0Xs′(Y,t)dY                =∫0X[1+ϵ0cos(2πL(Y−ct))]dY                =X+∫−2πctL2πX−ctLL2πϵ0cos(y)dy                =X+ϵ0L2π[sin(2π(X−ctL))+sin(2πctL)]

and, by differentiating with respect to time,

s˙(X,t)=cϵ0[cos(2πctL)−cos(2π(X−ct)L)].

Finally, in view of (5),

x˙1(t)=−ϵ0ccos(2πctL)                +c∫0Lϵ(X,t)(1+ϵ(X,t))1−pdX∫0L(1+ϵ(X,t))1−pdX

whence

x1(t)=−ϵ0L2πsin(2πLct)                +c∫0t∫0Lϵ(X,z)(1+ϵ(X,z))1−pdX∫0L(1+ϵ(X,z))1−pdXdz.

The Newtonian case is recovered by setting p = 0, i.e.,

x1(t)=−ϵ0L2πsin(2πLct)+ϵ02c2t.

Figure 4 displays three numerical examples. For p < 1 and, in particular for p = 0, the case of Newtonian resistance, we always have prograde motion (i.e., motion in the same direction as the one of the waves). This is indeed observed for example in snails, although in this case the force-velocity laws that we use in this paper would not be fully adequate to capture the properties of the mucus present between the animal and the surface [non-Newtonian rheology, suction effects, see (Denny, 1980) and (DeSimone et al., 2013)]. For p>1 and, in particular, for the limit case p = ∞ describing the perfect-grip/free-slip ideal version of the modulated friction laws typical of animals with setae, the motion is retrograde (i.e., motion in the opposite direction as the one of the waves). This is the behavior typically observed for earthworms.

Plot of x1(t) for a smooth contraction wave (26) for selected values of parameter p. The other parameters are ϵ0 = 0.6, L = 1 and c = 1.5.

 Square contraction wave

Consider the square contraction wave

                  ϵ(X,t):=ϵ0(X−ct)where ϵ0(x):={δif x~L≤ξ−δif x~L>ξ

or equivalently, in terms of the stretch,

s′(X,t)=1+ϵ(X,t)                 ={1+δif (X−ct)~L≤ξ1−δif(X−ct)~L>ξ

where L is the reference length of the crawler, c is the wave speed, ξ is the measure of the interval where ϵ = δ and the subscript ~L denotes the “modulo L” operator (i.e., y~L stands for ymodL).

By integrating the stretch over space, we get

s(X,t)=s(0,t)+∫0Xs′(Y,t)dY               =X+∫0Xϵ0(Y−ct)dY

whence

s˙(X,t)={−2δcif ct~L≤L−ξ ∧ X∈[ct,ct+ξ]~L2δcif ct~L>L−ξ ∧ X∈[ct+ξ−L,ct]~L0else

Finally, in view of (5),

x˙1(t)={A(p)if ct~L≤L−ξB(p)otherwise

where

A(p):=2δc(1+δ)1−pξ(1+δ)1−pξ+(1−δ)1−p(L−ξ)

and

B(p):=2δc(1−δ)1−p(ξ−L)(1+δ)1−pξ+(1−δ)1−p(L−ξ).

Defining α:=L-ξc and β:=Lc, it follows that

x1(t)=x1(0)+∫0tx1.(z)dz           =∫0⌊tβ⌋βx1.(z)dz+∫⌊tβ⌋β(tβ)βx1.(z)dz           =⌊tβ⌋[αA(p)+(β−α)B(p)]+           {{tβ}βA(p)if {tβ}β≤ααA(p)+({tβ}β−α)B(p)else

where {·} and ⌊·⌋ denote the fractional part and floor function respectively. The Newtonian case can be obtained as particular case by setting p = 0. Figure 5 shows three numerical examples and, as for the smooth contraction wave, the motion is prograde or retrograde whether p < 1 or p>1, respectively.

Plot of x1(t) for a square contraction wave (27) for selected values of parameter p. The other parameters are δ = 0.6, L = 1, T = 0.5 and c = 1.5.

 3.2. Discrete model: peristalsis as optimal gait

In the discrete framework, peristalsis is the result of phase coordination among the harmonic contractions of body segments, i.e., it has the form

ϵn(t)=ϱsin(2πtT+nΔφ) forn=1,…,N

where T is the period, ϱ is the amplitude and Δφ is the constant phase difference. As for the continuous case, discrete peristalsis produces prograde or retrograde motions according to the value of the parameter p in (4).

In this section we work out explicitly the problem of maximizing the displacement for a particular case from which peristalsis emerges, modulo an edge-effect.

Dissipation energy

Let us define an energy functional D:C2(ℝ, ℝN) → ℝ as

D[ϵ]:=∫0Td(t,ϵ,ϵ˙)dt

where

             d(t,ϵ,ϵ˙):=d1(t,ϵ,ϵ˙)+wd2(t,ϵ˙),d1(t,ϵ,ϵ˙):=∫0NL−1μfref(X,t)v(χ(X,t))dX,                 d2(t,ϵ˙):=∑n=1Nϵ˙n2(t)

i.e., the energy cost is the time integral over a period of a dissipation rate which is sum of two terms: d1(t,ϵ,ϵ°) is 1μ times the energy expended to overcome the friction force and d2(t,ϵ°) is the cost of control weighted by a scalar factor w. D[ϵ] is thus 1μ times the sum of the work due to the friction force plus the L2-norm of the controls suitably weighted to time the input direction.

Some calculations [see section 1 in Appendix C (Supplementary Material)] lead to

d1(t,ϵ,ϵ˙)=ϵ˙·D(ϵ)ϵ˙

where 𝔻(ϵ)∈ℝN×N for any ϵ∈(−1, +∞]N, and

d2(t,ϵ˙)=ϵ˙·INϵ˙

where 𝕀N is the N-dimensional identity matrix. Therefore the energy functional is

D[ϵ]=∫0Tϵ˙·G(ϵ)ϵ˙dt

where 𝔾(ϵ): = 𝔻(ϵ)+wIN.

 Non-linear optimal control problem

The non-linear optimization problem associated with energy functional (28) is

            maxϵ∈S Δx0:=∫0Tv(ϵ)·ϵ˙dtS={ϵ∈C2 | ϵ(0)=ϵ(T) ∧ D[ϵ]=c}.

The Euler-Lagrange equations lead to a second order non-linear system of ODEs, i.e., for n = 1, …, N,

ddt∂ℱ∂ϵ˙n(t,ϵ,ϵ˙)−∂ℱ∂ϵn(t,ϵ,ϵ˙)=0

where F(t,ϵ,ϵ°):=v(ϵ)·ϵ°-λϵ°·G(ϵ)ϵ°, λ is the Lagrange multiplier.

 The small-deformation regime

In the regime of small deformations we can expand the terms of problem (29) at the leading orders about ϵ = 0. As before, the net displacement per time period can be approximated by

Δx0≃V[ϵ,ϵ˙]:=∫0TVϵ·ϵ˙dt.

and the energy functional by

D[ϵ]=∫0Tϵ˙·G(ϵ)ϵ˙dt≃∫0Tϵ˙·Gϵ˙dt

where 𝔾: = 𝔾(0). Hence, in the small-deformation regime, the problem fits the form (17)-(18) for 𝔸 = 0 and 𝔹 = 𝔾. Moreover, G is bisymmetric (namely, symmetric about both of its diagonals) and depends only on N, L and w. Indeed

{G}ij={L34N(2i−1)(2(N−j)+1)if i<jL312N[4N(3i−2)−3(2i−1)2]+wif i=jL34N(2j−1)(2(N−i)+1)if i>j

[see section 2 in Appendix C (Supplementary Material)]. Therefore a solution must be of the form (25)

ϵn⋆(t)=ϱaϱnsin(2πTt+ϑn).

The centrosymmetry of 𝔾 and the skew-centrosymmetry of 𝕍 imply a reflectional symmetry about the center [see section 3 in Appendix B (Supplementary Material)], i.e.,

the moduli of components of e are symmetric about the center (cf. Figure 6), i.e., ϱN+1−n=ϱn ∀ n=1,…,N;

phase differences between adjacent segments are symmetric about the center, i.e., ϑn+1−ϑn=ϑN+1−n−ϑN−n ∀ n=1,…,N

so that the N-th phase differs from the (N−1)-th one by the same amount by which the second phase differs from the first one and so on (cf. Figure 6).

Plot of arguments and moduli of ϵn for n = 1, …, 15: amplitudes, (A), approximation by a IPD (Identical Phase Difference) model, (B), and relative errors, (C). Parameters: p = 100, w = 1, T = 1 and L = 1.

Equation (25) shows that the optimal gait requires a precise “phase coordination” of locomotion patterns among the segments, which is a common observation in Biology for several kinds of animals.

Numerical simulations show that the optimal solution turns out to be a discrete approximation of a traveling wave. In particular,

the moduli of en for n = 1, …, N can be approximated by a constant average value (cf. Figure 6), i.e., ϱn≃ϱ¯ constant

so that each segment undergoes a harmonic deformation with a certain initial phase;

phase differences between adjacent segments turn out to be almost constant, i.e., for a suitable φ0 ϑn≃nϑ⋆+ϑ0

holds true for n = 1, …, N (cf. Figure 6).

Therefore, in view of (32) and (33), the solution is a discrete approximation of a continuous traveling wave, i.e.,

ϵn⋆(t)=ϵ⋆(Yn,t)

where Yn=Xn+Xn+12 is the midpoint of the n-th segment and

ϵ⋆(X,t)=ϱaϱ(X)sin(2πTt+ϑ(X))                  ≃ϱaϱ¯sin(2πTt+ϑ⋆X+ϑ0)                  =ϱaϱ¯sin(ϑ⋆[X−(−1ϑ⋆2πT)t]+ϑ0)

whence the canonical form of traveling wave which describes peristalsis (cf. Figure 7)

ϵ⋆(X,t)≃H(X−vt)

where H(y):=ϱaϱ̄sin(ϑ⋆y+ϑ0) and v:=-1ϑ⋆2πT.

Plot of piecewise constant optimal strain ϵn⋆(X,t): the value is determined by means of the color legend. Parameters: p = 100, w = 10, T = 1, L = 1; N = 25.

 The edge-effect

The symmetric structure of the optimal gait (in the small-deformation regime) arises from underlying physical symmetries which clearly stand out in the properties of the matrices 𝔾 and 𝕍. In particular, an “edge-effect” is apparent: the 1D crawler is symmetric about its geometric center and segments near the edges behave differently with respect to adjacent segments, but in the same way as their centrosymmetric counterparts.

As expected, this edge-effect vanishes when considering an “infinite” (periodic) 1D crawler because, due to the shift-invariance symmetry, each segment behaves as a “geometric center.”

To show this claim consider a 1D crawler made up of infinitely many segments and assume that it is a periodic structure of which each module consists of N components (cf. Figure 8).

Kinematics of a discrete infinite 1D crawler consisting of identical segments of reference length L. (A) Reference configuration. (B) Current configuration.

At any time t, we already defined the relative displacement u(·, t) as the change of position of the material point X in the body's reference, i.e.,

x(t)=χ(X,t)=x0(t)+X+u(X,t).

The hypothesis of periodicity leads to

u(X+NL,t)=u(X,t) ∀ X,t.

From (35) we obtain that the friction force is periodic and we can consider the force balance in a single module. Therefore, condition (35) reads

∫0NLϵ(X,t)dX=0 ∀ t

and, in the discrete framework (9), this leads to

∑n=1Nϵn(t)=0 ∀ t.

The optimal control problem becomes

                 maxϵ∈S′ V[ϵ,ϵ˙]:=∫0Tϵ˙·VϵdtS′={ϵ∫0T∈C2(ℝ,ℝN) | ∑n=1Nϵn=0 ∧                    ϵ(0)=ϵ(T) ∧ ∫0Tϵ˙·Gϵ˙dt=c}

and it can be proved [see section 3 in Appendix C (Supplementary Material)] that its solutions need to be like (22), where the complex N-dimensional vector e has the form

e=[e1e2⋮en⋮eN]=[e1ei2πkNe1⋮ei2πkN(n−1)e1⋮ei2πkN(N−1)e1]

for some k∈{1, …, N−1} and e1∈ℂ\{0}, cf. Figure 9. In particular,

- each component of e has modulus ϱ: = ||e1||;

- each component can be obtained from the previous one by a rotation of 2πkN or, in other words, the phase difference between two consecutive components is constant, i.e., for n = 1, …, N arg(en)=(n−1)2πkN+arg(e1)=nϑ⋆+ϑ0

where ϑ⋆:=2πkN and ϑ0:=arg(e1)-2πkN;

Complex components of the vector e in the general case, (A), and in the periodic one, (B), for N = 8 segments. Parameters: p = 100, w = 10, T = 1, L = 1.

whence the exact harmonic peristalsis.

Notice that problem (38) can be written in terms of relative displacements un through the periodic version of transformation (15), i.e.,

ϵ(t)=Jperu(t)

where

Jper:=1L[1  −1−11   ⋱⋱   −11].

In particular, we get

                  maxϵ∈Su⋆ V[u,u˙]:=∫0Tu˙·Vu⋆udtSu⋆={u∈∫0TC3(ℝ,ℝN) | u(0)=u(T)                          ∧  E[u,u˙]:=∫0Tu˙·Gu⋆u˙dt=c}

where

Vu⋆:=JperTVJperGu⋆:=JperTGJper

are circulant matrices (namely, Toeplitz matrices where each row vector is rotated one element to the right relative to the preceding row vector), thus reflecting the geometric symmetry of the periodic structure, namely, the shift-invariance.

Considering the general (i.e., non-periodic) problem in terms of relative displacements yields

                       maxϵ∈Su V[u,u˙]:=∫0Tu˙·VuudtSu={∫0Tu∈C3(ℝ,ℝN)∫0T| u(0)=u(T)                               ∧  E[u,u˙]:=∫0Tu˙·Guu˙dt=c}

where

Vu:=JTVJGu:=JTGJ

are “quasi-circulant” matrices indeed

Vu:=Vu⋆+EVGu:=Gu⋆+EG

where 𝔼V and 𝔼G are null a part from the last column and the last row, i.e.,

and

where a = 2L3(3−N)+12Nw and b = L3(9−4N)−12Nw.

 Wavenumber

In the “periodic case” we can study the wavenumber (that is the number of waves travelling along the body of the crawler) of the optimal gait in relation to the number of metameres N and to the weight w. We fix the dissipation

E[ϵ,ϵ˙]:=∫0T(d1(t,ϵ,ϵ˙)+wd2(t,ϵ,ϵ˙))dt=c¯

and we let w vary from 0 to 102 for N∈[3, 250] (see Figure 10). As shown in Section 3 in Appendix C (Supplementary Material), the wavenumber of the optimal gait must be an integer close to the real number N2πarccos(126w-L33w+L3) and hence, for any fixed N, it depends on the order of magnitude of the weight w and

for w → ∞, it tends to 1, corresponding to a single wave spanning the whole length NL;

for w = 0, it is close to N3, i.e., one full wave-length every three segments.

Wavenumber of optimal gaits as a function of N and w. The axis of w∈[0, 100] is plotted on a log-scale with base 10. The color-bar gives the wavenumber. L = 1 and N∈[3, 250].

This behavior is qualitatively unaffected by the type of friction model which is adopted (i.e., by the choice of the parameter p).

 4. Discussion 4.1. Comparison with previous studies

To put our study in perspective, we consider the discrete framework and we compare our results with the ones presented by Fang et al. (2015). Here the authors perform an optimization of the so-called “average steady-state velocity” us among harmonic shape functions having the form (in our notation)

ϵn(t)=asin(2πTt+ηn) for n=1,…,N

where a∈(0,1L) is the oscillation amplitude, T is the period and ηn is the actuation phase for the n-th segment (or actuator).

Since the average steady-state velocity is given by

us=Δx0T=1T∫0Tv(ϵ)·ϵ˙dt ,

the optimization problem reads

maxη∈[0,2π)N us(η)=1T∫0Tv(ϵ)·ϵ˙dt

and in the small-deformation regime it can be replaced by

maxη∈[0,2π)N ∫0Tϵ˙·Vϵ dt.

Denote the actuation phase differences between adjacent segments by

pn:=Δη(n)=ηn+1−ηn for n=1,…,N−1.

From observations of numerical simulations, Fang et al. (2015) report that “[…] the optimized phase-different patterns are always reflectionally symmetric [about the center, Ed.] regardless of the initial symmetry requirements […]” and of the number of segments. Thus, a solution to (45) fulfills

pn=pN−n ∀ n.

In fact these properties can be rigorously proved under the assumption that problem (45) admits a unique solution in [0, 2π)N (see Appendix D in Supplementary Material).

Furthermore, property (47) can be proved also for (46), assuming it admits a unique solution in [0, 2π)N. To this aim, denote the unique solution to (46) by

ϵ˜(t):={ϵ˜n(t):=asin(2πTt+η˜n)}n=1,…,N,

and consider the shape change ϵ^(t) associated with

η^:=−Kη˜+2π

where

K:=[00⋯0100⋯10⋮⋮⋱⋮⋮01⋯0010⋯00]∈ℝN×N .

Notice that for n = 1, …, N,

ϵ^n(t):=asin(2πTt−(Kη˜)n)={−Kϵ˜(−t)}n

and hence, by exploiting the fact that V is skew-centrosymmetric (namely, 𝕂T𝕍𝕂 = −𝕍),

∫0Tϵ^.·Vϵ^dt=−∫0Tϵ˜.(−t)·KTVKϵ˜(−t)dt                      =∫−T0ϵ˜.·Vϵ˜dt=∫0Tϵ˜.·Vϵ˜dt.

Thus

η˜=−Kη˜+2π

which leads to (47).

Problem (46) constrains the L2-norm of the time-derivatives, i.e., for strains having the form (43) we get

∫0Tϵ˙·ϵ˙dt=2NT(aπ)2=:c⋆

regardless of η. Therefore we can extend the maximization to the C2 periodic strains whose time derivative fulfills the same constraint, i.e.,

                       maxϵ∈S V[ϵ,ϵ˙]:=∫0Tϵ˙·Vϵ dtS={∫0Tϵ∈C2(ℝ,ℝN) | ϵ(0)=ϵ(T)∧ ∫0Tϵ˙·ϵ˙dt=c⋆}.

Since problem (21) reduces to (48) when 𝔸 = 0 and 𝔹 = 𝕀N, a solution to (48) must be of the form

ϵn⋆(t)=aN ||en|| sin(2πTt+arg(en)+ϑa)

where e = (en)n is a unit eigenvector associated with the maximum-modulus eigenvalue of V and ϑa is a constant. Notice that the reflectional symmetry about the center still holds true. As a matter of fact, (49) leads to a slight increment in the net displacement with respect to the solution to (46), cf. Figure 11.

Average velocities us, (44) obtained by the solution ϵ~ to (46) (blue bars) and by the solution ϵ⋆ to (48) (yellow bars) for different numbers of segments: N = 25, 50, 100, 150, 200. (A) is for p = 0 (Newtonian case) and (B) for p = 100. The other parameters are T = 2π, a = 2−10, L = 1.

 4.2. Summary and outlook

Our analysis confirms the effectiveness of mimicking peristalsis in bio-inspired robots, at least in the small-deformation regime. This bio-inspired actuation strategy has been implemented on a trial-and-error basis many times in the robotics literature and, more recently, also proposed as optimal (in some suitably defined sense, and in some suitably defined class of actuation strategies). Our main result is a mathematically rigorous proof that, in the small deformation regime, actuation by peristaltic waves is an optimal control strategy emerging naturally from the geometric symmetry of the system, namely, the invariance under shifts along the body axis. This is true exactly in the periodic case, and approximately true in the case of finite length, modulo edge-effects. Our results is of theoretical nature. Nevertheless, we believe that it has important consequences for applications. For example, it helps us not to fall into the naive temptation to expect that peristaltic waves are always an optimal actuation strategy just because they are observed in nature, but to exercise critical judgment whenever the hypotheses on the geometric symmetry that are “responsible” for the optimality result in our case (invariance under shift of a homogeneous one-dimensional system) are false.

Actuation by phase coordination, optimal actuation by identical phase difference, and the connections between this and traveling waves have been already discussed in the literature (e.g., Fang et al., 2015), but never through a mathematically rigorous analysis of the optimal control problem, of the symmetry properties of the governing equations and operators, and of the relation between these and the geometric symmetries of the system. This is exactly what we do in this paper. The added value of this analysis is that we are able to show (for the first time, to the best of our knowledge, at least in the robotics literature) that peristaltic waves are the signature of the invariance with respect to shifts (a geometric symmetry) of a homogeneous one-dimensional system.

Further work will be needed to test the effectiveness of peristaltic waves as a locomotion strategies if large deformations are allowed. In addition, future work will explore the issue of how peristalsis is actually enforced in biological systems. Of particular interest is the dichotomy between the paradigm of actuations via a Central Pattern Generator (CPG), as opposed to local sensory and feedback mechanisms. The CPG paradigm is apparent in several different organisms (Marder and Bucher, 2001; Grillner, 2006) and has been employed in robotics with some success (Ijspeert, 2008; Boxerbaum et al., 2012). However, there is a growing awareness of the role played by proprioception, especially for lower organisms such as the nematode worm C. elegans (Boyle et al., 2012; Wen et al., 2012) and D. melanogaster larvae (Pehlevan et al., 2016).

 Author contributions

AD and FA conceived research. DA executed research and performed numerical simulation. AD supervised research. FA contributed expertise on circulant matrices. All authors analyzed the data and wrote the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The handling Editor declared a shared affiliation, though no other collaboration, with one of the authors, AD.

This article was written while DA and AD where visiting Division C of the Engineering Department of the University of Cambridge. We thank Prof. Vikram Deshpande and Corpus Christi College for their warm hospitality.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frobt.2018.00099/full#supplementary-material

 References Alberts B. Johnson A. Lewis J. Raff M. (2002). Molecular Biology of the Cell, 4th Edn. New York, NY: Garland Science. Boxerbaum A. Shaw K. Chiel H. Quinn R. (2012). Continuous wave peristaltic motion in a robot. Int. J. Robot. Res. 31, 302–318. 10.1177/0278364911432486 Boyle J. H. Berri S. Cohen N. (2012). Gait modulation in c. elegans: an integrated neuromechanical model. Front. Comput. Neurosci. 6:10. 10.3389/fncom.2012.0001022408616 Collar A. (1962). On centrosymmetric and centroskew matrices. Q. J. Mechan. Appl. Math. 15, 265–281. 10.1093/qjmam/15.3.265 Daltorio K. Boxerbaum A. Horchler A. S. Shaw K. M. Chiel H. J. Quinn R. D. (2013). Efficient worm-like locomotion: slip and control of soft-bodied peristaltic robots. Bioinspir. Biomimet. 8:035003. 10.1088/1748-3182/8/3/03500323981561 Denny M. (1980). The role of gastropod pedal mucus in locomotion. Nature 285:160. 10.1038/285160a0 DeSimone A. Guarnieri F. Noselli G. Tatone A. (2013). Crawlers in viscous environments: linear vs non-linear rheology. Int. J. Non Linear Mech. 56, 142–147. 10.1016/j.ijnonlinmec.2013.02.007 DeSimone A. Tatone A. (2012). Crawling motility through the analysis of model locomotors: two case studies. Eur. Phys. J. E Soft Matt. Biol. Phys. 35, 1–8. 10.1140/epje/i2012-12085-x22972227 Edwards C. Hendrix P. Arancon N. (2018). Biology and Ecology of Earthworms. Springer. Available online at: https://www.springer.com/in/book/9780387749426 Fang H. Li S. Wang K. W. Xu J. (2015). Phase coordination and phase–velocity relationship in metameric robot locomotion. Bioinspir. Biomimet. 10:066006. 10.1088/1748-3190/10/6/066006 Gardner C. (1976). The neuronal control of locomotion in the earthworm. Biol. Rev. 51, 25–52. 10.1111/j.1469-185X.1976.tb01119.x766843 Garrey W. Moore A. (1915). Peristalsis and coordination in the earthworm. Am. J. Physiol. Legacy Cont. 39, 139–148. 10.1152/ajplegacy.1915.39.2.139 Ge J. Calderón A. Pérez-Arancibia N. (2017). An earthworm-inspired soft crawling robot controlled by friction. arXiv [preprint] arXiv:1707.04084. Gray J. Lissmann H. (1938). Studies in animal locomotion: Vii. locomotory reflexes in the earthworm. J. Exp. Biol. 15, 506–517. Grillner S. (2006). Biological pattern generation: the cellular and computational logic of networks in motion. Neuron 52, 751–766. 10.1016/j.neuron.2006.11.00817145498 Ijspeert A. (2008). Central pattern generators for locomotion control in animals and robots: a review. Neural Netw. 21, 642–653. 10.1016/j.neunet.2008.03.01418555958 Marder E. Bucher D. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, R986–R996. 10.1016/S0960-9822(01)00581-411728329 McInerney A. (2013). First Steps in Differential Geometry: Riemannian, Contact, Symplectic. Undergraduate Texts in Mathematics. New York, NY: Springer. Menciassi A. Accoto D. Gorini S. Dario P. (2006). Development of a biomimetic miniature robotic crawler. Auton. Robots 21, 155–163. 10.1007/s10514-006-7846-9 Nemitz M. P. Mihaylov P. Barraclough T. W. Ross D. Stokes A. A. (2016). Using voice coils to actuate modular soft robots: wormbot, an example. Soft Robot. 3, 198–204. 10.1089/soro.2016.000928078195 Pehlevan C. Paoletti P. Mahadevan L. (2016). Integrative neuromechanics of crawling in D. melanogaster larvae. Elife 5:23. 10.7554/eLife.11031 Quillin K. (1999). Kinematic scaling of locomotion by hydrostatic animals: ontogeny of peristaltic crawling by the earthworm lumbricus terrestris. J. Exp. Biol. 202, 661–674. 10021320 Umedachi T. Kano T. Ishiguro A. Trimmer B. (2016). Gait control in a soft robot by sensing interactions with the environment using self-deformation. R. Soc. Open Sci. 3:160766. 10.1098/rsos.16076628083114 Van Brunt B. (2004). The Calculus of Variations. New York, NY: Springer. Wang K. Yan G. Ma G. Ye D. (2009). An earthworm-like robotic endoscope system for human intestine: design, analysis, and experiment. Anna. Biomed. Eng. 37, 210–221. 10.1007/s10439-008-9597-619003537 Wen Q. Po M. D. Hulme E. Chen S. Liu X. Kwok S. W. . (2012). Proprioceptive coupling within motor neurons drives c. elegans forward locomotion. Neuron 76, 750–761. 10.1016/j.neuron.2012.08.03923177960 Wiezel O. Giraldi L. DeSimone A. Or Y. Alouges F. (2018). Energy-optimal small-amplitude strokes for multi-link microswimmers: Purcell's loops and taylor's waves reconciled. arXiv [preprint] arXiv:1801.04687.

Funding. We gratefully acknowledge the financial support of the European Research Council through the Advanced Grant 340685-MicroMotility.