Peristaltic Waves as Optimal Gaits in Metameric Bio-Inspired Robots

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.

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

$\left[X_1:=0,X_2:=L\right],$

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)=\chi (X,t):=x_1(t)+s(X,t)\in \left[x_1(t),x_2(t)\right]$

where x₁(t): = χ(X₁, t) and x₂(t): = χ(X₂, t). By definition, s(0, t)≡0 for all t and we assume that

$\begin{array}{cc}{s}^{\prime }(X,t):=\frac{\partial s(X,t)}{\partial X}>0\forall X,t& (1)\end{array}$

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

FIGURE 1

Figure 1. 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

$u_a(X,t):=\chi (X,t)-X$

so that

${\chi }^{\prime }(X,t):=\frac{\partial \chi (X,t)}{\partial X}=1+u_{a^{\prime }}(X,t)=1+ϵ(X,t)$

where we have implicitly defined the strain

$ϵ(X,t):=\frac{\partial u_a}{\partial X},$

in terms of which condition (1) reads

$\begin{array}{cc}ϵ(X,t)>-1\forall X,t.& (2)\end{array}$

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

$\dot{\chi }(X,t)={\dot{x}}_1(t)+\dot{s}(X,t),$

the spatial (or Eulerian) velocity is given by

$v(x,t)=\dot{\chi }(X,t){|}_{X={\chi }^{-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 x₁(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

$\begin{array}{cc}f(x,t):=-\mu v(x,t)& (3)\end{array}$

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

• or a more general “p-model”

$\begin{array}{l}f(x,t):=-\mu g_p(ϵ(X,t))v(x,t)\\ \text{ where\hspace{1em}}{g}_p(ϵ):={\left(\frac{1}{1+ϵ}\right)}^p& (4)\end{array}$

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 g_p(ϵ) 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).

FIGURE 2

Figure 2. Function g_p(ϵ) 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.,

$F_f(t)={\displaystyle {\int }_{x_1}^{x_2}f}(x,t)dx={\displaystyle {\int }_0^L{f}_{ref}(X,t)}dX$

where

$f_{ref}(X,t):=f\left(x_1(t)+s(X,t),t\right)s^{\prime }(X,t).$

Then, by neglecting inertia, the force balance yields

$\begin{array}{l}0=F_f(t)+F_e(t)\\ =-\mu {\displaystyle {\int }_0^L{\left(\frac{1}{1+ϵ(X,t)}\right)}^p}\left({\dot{x}}_1(t)+\dot{s}(X,t)\right)s^{\prime }(X,t)dX+F_e(t)\\ =\left[-\mu {\displaystyle {\int }_0^L{\left(\frac{1}{1+ϵ(X,t)}\right)}^p}{s}^{\prime }(X,t)dX\right]{\dot{x}}_1\\ -\mu {\displaystyle {\int }_0^L{\left(\frac{1}{1+ϵ(X,t)}\right)}^p}\dot{s}(X,t)s^{\prime }(X,t)dX+F_e(t)\\ =\left[-\mu {\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}dX\right]{\dot{x}}_1\\ -\mu {\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}\dot{s}(X,t)dX+F_e(t).\end{array}$

The square bracket multiplying ẋ₁(t) in the formula above is the drag for rigid motion at unit speed and fixed shape ϵ(X, t), while F_e(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 ẋ₁(t), we obtain

$\begin{array}{l}{\dot{x}}_1(t)=-\frac{{\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}\dot{s}(X,t)dX}{{\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}dX}\\ +\frac{F_e(t)}{\mu {\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}dX}& (5)\end{array}$

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

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

$\begin{array}{l}s(X,t)=s(0,t)+{\displaystyle {\int }_0^X{s}^{\prime }}(Y,t)dY\\ ={\displaystyle {\int }_0^Xϵ}(Y,t)dY+X& (6)\end{array}$

one has

$\begin{array}{cc}\dot{s}(X,t)={\displaystyle {\int }_0^X\dot{ϵ}}(Y,t)dY& (7)\end{array}$

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)

$\left(X_{n-1}:=(n-1)L,X_n:=nL\right)$

for n = 1, …, N. Let x₀(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):=\chi (X,t)=x_0(t)+s(X,t)=X+u_a(X,t)$

whence the definition of strain

$\begin{array}{cc}ϵ(X,t):=u_{a^{\prime }}(X,t)=s^{\prime }(X,t)-1.& (8)\end{array}$

FIGURE 3

Figure 3. 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.,

$\begin{array}{cc}ϵ(X,t):=\{\begin{array}{ll}\begin{array}{l}{ϵ}_1(t)\\ \end{array}\hfill & \begin{array}{l}\text{if }X\in \left[X_0,X_1\right)\\ (\text{segment 1})\end{array}\hfill \\ ⋮\hfill & \hfill \\ \begin{array}{l}{ϵ}_N(t)\\ \end{array}\hfill & \begin{array}{l}\text{if }X\in \left(X_{N-1},X_N\right]\\ (\text{segment N})\end{array}\hfill \end{array}& (9)\end{array}$

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

$\begin{array}{cc}\dot{s}(X,t)=L{\displaystyle \sum _{j=1}^{n-1}{\dot{ϵ}}_j}(t)+\left[X-(n-1)L\right]{\dot{ϵ}}_n(t)& (10)\end{array}$

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

${ϵ}_n(t)>-1\text{\hspace{1em}for all }t\text{ and for }n=1,\dots ,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

$\begin{array}{l}0=\left[-\mu {\displaystyle {\int }_0^{NL}{\left(1+ϵ(X,t)\right)}^{1-p}}dX\right]{\dot{x}}_0(t)\\ -\mu {\displaystyle {\int }_0^{NL}{\left(1+ϵ(X,t)\right)}^{1-p}}\dot{s}(X,t)dX+F_e(t)\end{array}$

where, in view of (9),

${\displaystyle {\int }_0^{NL}{\left(1+ϵ(X,t)\right)}^{1-p}}dX=L{\displaystyle \sum _{j=1}^N{(1+{ϵ}_j)}^{1-p}}$

and, in view of (10),

$\begin{array}{l}{\displaystyle {\int }_0^{NL}{\left(1+ϵ(X,t)\right)}^{1-p}}\dot{s}(X,t)dX\\ ={\displaystyle \sum _{n=1}^N{(1+{ϵ}_n)}^{1-p}}·{\displaystyle {\int }_{(n-1)L}^{nL}\left[\text{}\left(X-(n-1)L\right){\dot{ϵ}}_n(t)+L{\displaystyle \sum _{k=1}^{n-1}{\dot{ϵ}}_k}(t)\text{}\right]}dX\\ =\frac{L^2}{2}{\displaystyle \sum _{n=1}^N{(1+{ϵ}_n)}^{1-p}}\left[{\dot{ϵ}}_n(t)+2{\displaystyle \sum _{k=1}^{n-1}{\dot{ϵ}}_k}(t)\right]\\ =\frac{L^2}{2}{\displaystyle \sum _{j=1}^N\left[{(1+{ϵ}_j)}^{1-p}+2{\displaystyle \sum _{m=j+1}^N{(1+{ϵ}_m)}^{1-p}}\right]}{\dot{ϵ}}_j(t).\end{array}$

Solving for ẋ₀(t) we obtain

$\begin{array}{l}{\dot{x}}_0(t)={\left[{\displaystyle \sum _{j=1}^N{(1+{ϵ}_j)}^{1-p}}\right]}^{-1}\{\text{}\frac{F_e(t)}{\mu L}\\ -\frac{L}{2}{\displaystyle \sum _{j=1}^N\left[\text{}{(1+{ϵ}_j)}^{1-p}+2{\displaystyle \sum _{m=j+1}^N{(1+{ϵ}_m)}^{1-p}}\text{}\right]}\text{}{\dot{ϵ}}_j(t)\text{}\}& (11)\end{array}$

which can be rewritten in the following vectorial form

$\begin{array}{l}{\dot{x}}_0(t)={\displaystyle \sum _{n=1}^N{v}_n}(ϵ(t)){\dot{ϵ}}_n(t)+q(ϵ(t))F_e(t)\\ =v(ϵ(t))·\stackrel{.}{ϵ}(t)+q(ϵ(t))F_e(t)\end{array}$

where

$\begin{array}{l}{v}_n(ϵ)=-\frac{L}{2}\frac{{\left(1+{ϵ}_n\right)}^{1-p}+2{\displaystyle \sum _{m=n+1}^N{\left(1+{ϵ}_m\right)}^{1-p}}}{{\displaystyle \sum _{j=1}^N{\left(1+{ϵ}_j\right)}^{1-p}}},\\ q(ϵ(t))={\left[\mu L{\displaystyle \sum _{j=1}^N{\left(1+{ϵ}_j(t)\right)}^{1-p}}\right]}^{-1},\\ v(ϵ)=\left[\begin{array}{c}{v}_1(ϵ)\\ ⋮\\ v_N(ϵ)\end{array}\right]\text{\hspace{1em}and\hspace{1em}}ϵ(t)=\left[\begin{array}{c}{ϵ}_1(t)\\ ⋮\\ {ϵ}_N(t)\end{array}\right].\end{array}$

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

$\begin{array}{l}\Delta x_0=x_0(T)-x_0(0)\\ ={\displaystyle {\int }_0^{T}v}(ϵ(t))·\stackrel{.}{ϵ}(t)dt+{\displaystyle {\int }_0^{T}q}(ϵ(t))F_e(t)dt.& (12)\end{array}$

Relative Displacement

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

$\begin{array}{l}u(X,t):=u_a(X,t)-x_0(t)\\ =s(X,t)-X\\ ={\displaystyle {\int }_0^Xϵ}(Y,t)dY& (13)\end{array}$

in order to describe the displacement in a coordinate system which is “co-moving” with the left end x₀(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∈[X_n−1, X_n],

$\begin{array}{cc}u(X,t)=L{\displaystyle \sum _{j=1}^{n-1}{ϵ}_j}(t)+[X-(n-1)L]{ϵ}_n(t).& (14)\end{array}$

Setting

$u_n(t):=u(X_n,t)\text{\hspace{1em}for }n=1,\dots ,N,$

we have

${ϵ}_n=\frac{u_n-u_{n-1}}{L}\text{\hspace{1em}for }n=1,\dots ,N,$

where u₀≡0. Equivalently

$\begin{array}{cc}ϵ(t)=\mathbb{J}u(t)& (15)\end{array}$

where $ϵ={({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_N)}^T$, $\text{u}={(u_1,u_2,\dots ,u_N)}^T$ and

$\begin{array}{cc}\mathbb{J}:=\frac{1}{L}\left[\begin{array}{cccc}1& & & \\ -1& 1& & \\ & \ddots & \ddots & \\ & & -1& 1\end{array}\right].& (16)\end{array}$

Optimal Control Problems

In this section we address the problem of maximizing the net displacement Δx₀ 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 F_e≡0.

Feasible Region

We assume that the shape function ϵ(t) is a C² 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\left[ϵ,\stackrel{°}{ϵ}\right]=c$, where the energy functional is assumed to be of the following quadratic form (in both ϵ and $\stackrel{°}{ϵ}$)

$\begin{array}{cc}E[ϵ,\stackrel{.}{ϵ}]:={\displaystyle {\int }_0^T\mathbb{A}}ϵ·ϵdt+{\displaystyle {\int }_0^T\mathbb{B}}\stackrel{.}{ϵ}·\stackrel{.}{ϵ}dt& (17)\end{array}$

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

$\mathcal{S}:=\{ϵ\in C^2\left(ℝ,{ℝ}^N\right)|ϵ(0)=ϵ(T)\wedge |E[ϵ,\stackrel{.}{ϵ}]=c\}.$

Optimization Problem

The general (non-linear) optimization problem is

$\begin{array}{cc}\underset{ϵ\in \mathcal{S}}{\max }\Delta x_0\left[ϵ\right]:={\displaystyle {\int }_0^{T}v}(ϵ(t))·\dot{ϵ}(t)dt& (18)\end{array}$

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,

$\begin{array}{cc}\frac{d}{dt}\frac{\partial ℱ}{\partial {\dot{ϵ}}_n}(t,ϵ,\dot{ϵ})-\frac{\partial ℱ}{\partial {ϵ}_n}(t,ϵ,\dot{ϵ})=0& (19)\end{array}$

where ${F}(t,ϵ,\stackrel{°}{ϵ}):=\text{v}(ϵ)·\stackrel{°}{ϵ}-\lambda (Aϵ·ϵ+B\stackrel{°}{ϵ}·\stackrel{°}{ϵ})$ (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.,

$\begin{array}{l}\Delta x_0={\displaystyle {\int }_0^{T}v(ϵ)·\dot{ϵ}}dt\\ ={\displaystyle {\int }_0^T\left(v(0)+v_{ϵ}(0)ϵ+o(ϵ)\right)·\dot{ϵ}}dt\\ \simeq v(0)·\left(ϵ(T)-ϵ(0)\right)+{\displaystyle {\int }_0^T{v}_{ϵ}(0)ϵ·\dot{ϵ}}dt\\ ={\displaystyle {\int }_0^T{v}_{ϵ}(0)ϵ·\dot{ϵ}}dt\end{array}$

and, integrating by parts,

$\begin{array}{l}{\displaystyle {\int }_0^T{v}_{ϵ}(0)ϵ·}\dot{ϵ}dt={\left[v_{ϵ}(0)ϵ·ϵ\right]}_0^T-{\displaystyle {\int }_0^T{v}_{ϵ}(0)\dot{ϵ}·ϵ}dt\\ ={\displaystyle {\int }_0^T-\dot{ϵ}·v_{ϵ}^T(0)ϵ}dt\end{array}$

whence

$\begin{array}{l}\Delta x_0\simeq \frac{1}{2}\left[{\displaystyle {\int }_0^T{v}_{ϵ}(0)ϵ·\dot{ϵ}}dt+{\displaystyle {\int }_0^T-v_{ϵ}^T(0)ϵ·\dot{ϵ}}dt\right]\\ ={\displaystyle {\int }_0^T\left(\frac{v_{ϵ}(0)-v_{ϵ}^T(0)}{2}\right)}ϵ·\dot{ϵ}dt=:V[ϵ,\dot{ϵ}].\end{array}$

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

${\left\{\mathbb{V}\right\}}_{ij}=\{\begin{array}{ll}L(p-1)\frac{i-j+N}{2N^2}\hfill & \text{if }i<j\hfill \\ 0\hfill & \text{if }i=j\hfill \\ -L(p-1)\frac{j-i+N}{2N^2}\hfill & \text{if }i>j\hfill \end{array}$

(see Appendix A in Supplementary Material).

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

$\begin{array}{cc}\underset{ϵ\in \mathcal{S}}{\max }V[ϵ,\dot{ϵ}]:={\displaystyle {\int }_0^T\dot{ϵ}}·\mathbb{V}ϵdt.& (20)\end{array}$

The corresponding Euler-Lagrange equations

$\frac{d}{dt}\frac{\partial ℒ}{\partial \dot{ϵ}}(t,ϵ,\dot{ϵ})-\frac{\partial ℒ}{\partial ϵ}(t,ϵ,\dot{ϵ})=0$

where

$ℒ(t,ϵ,\dot{ϵ}):=\mathbb{V}ϵ·\dot{ϵ}-\lambda \left(\mathbb{A}ϵ·ϵ+\mathbb{B}\dot{ϵ}·\dot{ϵ}\right),$

lead to the following system of second order linear ODEs

$\begin{array}{cc}\mathbb{V}\dot{ϵ}=\lambda \left(\mathbb{B}\ddot{ϵ}-\mathbb{A}ϵ\right).& (21)\end{array}$

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 ${𝔹}^{-\frac{1}{2}}𝕍{𝔹}^{-\frac{1}{2}}$ (resp. ${𝔸}^{-\frac{1}{2}}𝕍{𝔸}^{-\frac{1}{2}}$) have dimension 1. Indeed, as shown in sections 1 and 2 in Appendix B (Supplementary Material), it turns out that

• for 𝔸 = 0 and 𝔹 symmetric and positive definite, up to a constant, a solution of (20) must be of the form

$\begin{array}{cc}ϵ(t)=-\frac{T}{\pi }\Re \left(\alpha i{e}^{\frac{2\pi i}{T}t}e\right)& (22)\end{array}$

where α∈ℂ\{0} is a constant such that $\left|\right|\alpha \left|\right|=\sqrt{\frac{c}{2T}}$ and $\text{e}={(e_1,e_2,\dots ,e_N)}^T\in {ℂ}^N\backslash \left\{0\right\}$ is a suitable constant vector depending only on 𝔸 and 𝕍.

• for 𝔸 symmetric and positive definite and 𝔹 = 0, a solution of (20) with ϵ of unitary time frequency must be of the form

$\begin{array}{cc}ϵ(t)=2\Re \left(\alpha e^{\frac{2\pi i}{T}t}e\right)& (23)\end{array}$

where α∈ℂ\{0} is a constant such that $\left|\right|\alpha \left|\right|=\sqrt{\frac{c}{2T}}$ and $\text{e}={(e_1,e_2,\dots ,e_N)}^T\in {ℂ}^N\backslash \left\{0\right\}$ is a suitable constant vector depending only on 𝔸 and 𝕍.

Both (22) and (23) have the form

$ϵ(t)=\Re \left(\widehat{\alpha }{e}^{\frac{2\pi i}{T}t}e\right)$

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

$\widehat{\alpha }={\varrho }_a{e}^{i{\vartheta }_a}\text{\hspace{1em}and\hspace{1em}}{e}_n={\varrho }_n{e}^{i{\vartheta }_n}\forall n,$

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

$\begin{array}{l}{ϵ}_n(t)={\varrho }_a{\varrho }_n\Re \left(e^{i\left(\frac{2\pi }{T}t+{\vartheta }_a+{\vartheta }_n\right)}\right)\\ ={\varrho }_a{\varrho }_n\sin \left(\frac{2\pi }{T}t+{\vartheta }_a+{\vartheta }_n+\frac{\pi }{2}\right),& (24)\end{array}$

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 ${\vartheta }_a+\frac{\pi }{2}=0$, i.e.,

$\begin{array}{cc}{ϵ}_n(t)={\varrho }_a{\varrho }_n\sin \left(\frac{2\pi }{T}t+{\vartheta }_n\right),& (25)\end{array}$

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

$\begin{array}{l}\underset{ϵ\in \mathcal{S}}{\max }V[ϵ,\dot{ϵ}]\\ \mathcal{S}:=\left\{ϵ\in C^2|ϵ(0)=ϵ(T)\wedge {\displaystyle {\int }_0^T|}|ϵ|{|}_{{ℝ}^N}^{2}dt=1\right\}\end{array}$

then it is a solution also to

$\begin{array}{l}{\max }_{ϵ\in \mathcal{S}}V[ϵ,\dot{ϵ}]\\ \mathcal{S}:=\{ϵ\in C^2|ϵ(0)=ϵ(T){{\displaystyle \int }}_0^T\wedge {\displaystyle {\int }_0^T|}\left|\dot{ϵ}\right|{|}_{{ℝ}^N}^{2}dt={\left(\frac{2\pi }{T}\right)}^2\}\end{array}$

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

$\begin{array}{l}{\displaystyle {\int }_0^T\mathbb{A}}ϵ·ϵdt\ge {\lambda }_{min}(\mathbb{A}){\displaystyle {\int }_0^Tϵ}·ϵdt\\ ={\lambda }_{min}(\mathbb{A}){\left(\frac{T}{2\pi }\right)}^2{\displaystyle {\int }_0^T\dot{ϵ}}·\dot{ϵ}dt\\ \ge {\left(\frac{T}{2\pi }\right)}^2\frac{{\lambda }_{min}(\mathbb{A})}{{\lambda }_{max}(\mathbb{B})}{\displaystyle {\int }_0^T\mathbb{B}}\dot{ϵ}·\dot{ϵ}dt\end{array}$

and, analogously,

$\begin{array}{l}{\displaystyle {\int }_0^T\mathbb{B}}\dot{ϵ}·\dot{ϵ}dt\ge {\lambda }_{min}(\mathbb{B}){\displaystyle {\int }_0^T\dot{ϵ}}·\dot{ϵ}dt\\ ={\lambda }_{min}(\mathbb{B}){\left(\frac{2\pi }{T}\right)}^2{\displaystyle {\int }_0^Tϵ}·ϵdt\\ \ge {\left(\frac{2\pi }{T}\right)}^2\frac{{\lambda }_{min}(\mathbb{B})}{{\lambda }_{max}(\mathbb{A})}{\displaystyle {\int }_0^T\mathbb{A}}ϵ·ϵdt.\end{array}$

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., F_e≡0.

Smooth Contraction Wave

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

$\begin{array}{cc}ϵ(X,t):={ϵ}_0\cos \left(\frac{2\pi }{L}(X-ct)\right)& (26)\end{array}$

or equivalently, in terms of the stretch,

$s^{\prime }(X,t)=1+ϵ(X,t)=1+{ϵ}_0\cos \left(\frac{2\pi }{L}(X-ct)\right)$

where ϵ₀ 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,

$\begin{array}{l}s(X,t)=s(0,t)+{\displaystyle {\int }_0^X{s}^{\prime }}(Y,t)dY\\ ={\displaystyle {\int }_0^X\left[1+{ϵ}_0\cos \left(\frac{2\pi }{L}(Y-ct)\right)\right]}dY\\ =X+{\displaystyle {\int }_{-2\pi \frac{ct}{L}}^{2\pi \frac{X-ct}{L}}\frac{L}{2\pi }}{ϵ}_0\cos (y)dy\\ =X+\frac{{ϵ}_{0}L}{2\pi }\left[\sin \left(2\pi \left(\frac{X-ct}{L}\right)\right)+\sin \left(2\pi \frac{ct}{L}\right)\right]\end{array}$

and, by differentiating with respect to time,

$\dot{s}(X,t)=c{ϵ}_0\left[\cos \left(2\pi \frac{ct}{L}\right)-\cos \left(2\pi \frac{(X-ct)}{L}\right)\right].$

Finally, in view of (5),

$\begin{array}{l}{\dot{x}}_1(t)=-{ϵ}_{0}c\cos \left(\frac{2\pi ct}{L}\right)\\ +c\frac{{\displaystyle {\int }_0^Lϵ}(X,t){\left(1+ϵ(X,t)\right)}^{1-p}dX}{{\displaystyle {\int }_0^L{\left(1+ϵ(X,t)\right)}^{1-p}}dX}\end{array}$

whence

$\begin{array}{l}{x}_1(t)=-\frac{{ϵ}_{0}L}{2\pi }\sin \left(\frac{2\pi }{L}ct\right)\\ +c{\displaystyle {\int }_0^t\frac{{\displaystyle {\int }_0^Lϵ}(X,z){\left(1+ϵ(X,z)\right)}^{1-p}dX}{{\displaystyle {\int }_0^L{\left(1+ϵ(X,z)\right)}^{1-p}}dX}}dz.\end{array}$

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

$x_1(t)=-\frac{{ϵ}_{0}L}{2\pi }\sin \left(\frac{2\pi }{L}ct\right)+\frac{{ϵ}_0^{2}c}{2}t.$

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.

FIGURE 4

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

Square Contraction Wave

Consider the square contraction wave

$\begin{array}{l}ϵ(X,t):={ϵ}_0(X-ct)\\ \text{where\hspace{1em}}{ϵ}_0(x):=\{\begin{array}{l}\begin{array}{ll}\delta \hfill & \text{if }{x}_{~L}\le \xi \hfill \\ -\delta \hfill & \text{if }{x}_{~L}>\xi \hfill \end{array}\hfill \end{array}& (27)\end{array}$

or equivalently, in terms of the stretch,

$\begin{array}{l}{s}^{\prime }(X,t)=1+ϵ(X,t)\\ =\{\begin{array}{l}\begin{array}{ll}1+\delta \hfill & \text{if }{\left(X-ct\right)}_{~L}\le \xi \hfill \\ 1-\delta \hfill & \text{if}{\left(X-ct\right)}_{~L}>\xi \hfill \end{array}\hfill \end{array}\end{array}$

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

$\begin{array}{l}s(X,t)=s(0,t)+{\displaystyle {\int }_0^X{s}^{\prime }}(Y,t)dY\\ =X+{\displaystyle {\int }_0^X{ϵ}_0}(Y-ct)dY\end{array}$

whence

$\dot{s}(X,t)=\{\begin{array}{ll}-2\delta c\hfill & \text{if\hspace{1em}}c{t}_{~L}\le L-\xi \hfill \\ \hfill & \wedge X\in {\left[ct,ct+\xi \right]}_{~L}\hfill \\ 2\delta c\hfill & \text{if\hspace{1em}}c{t}_{~L}>L-\xi \hfill \\ \hfill & \wedge X\in {\left[ct+\xi -L,ct\right]}_{~L}\hfill \\ 0\hfill & \text{else}\hfill \end{array}$

Finally, in view of (5),

${\dot{x}}_1(t)=\{\begin{array}{ll}A(p)\hfill & \text{if }c{t}_{~L}\le L-\xi \hfill \\ B(p)\hfill & \text{otherwise}\hfill \end{array}$

where

$A(p):=\frac{2\delta c{(1+\delta )}^{1-p}\xi }{{(1+\delta )}^{1-p}\xi +{(1-\delta )}^{1-p}(L-\xi )}$

and

$B(p):=\frac{2\delta c{(1-\delta )}^{1-p}(\xi -L)}{{(1+\delta )}^{1-p}\xi +{(1-\delta )}^{1-p}(L-\xi )}.$

Defining $\alpha :=\frac{L-\xi }{c}$ and $\beta :=\frac{L}{c}$, it follows that

$\begin{array}{l}{x}_1(t)=x_1(0)+{\displaystyle {\int }_0^t\stackrel{.}{x_1}}(z)dz\\ ={\displaystyle {\int }_0^{\lfloor \frac{t}{\beta }\rfloor \beta }\stackrel{.}{x_1}}(z)dz+{\displaystyle {\int }_{\lfloor \frac{t}{\beta }\rfloor \beta }^{\left(\frac{t}{\beta }\right)\beta }\stackrel{.}{x_1}}(z)dz\\ =\lfloor \frac{t}{\beta }\rfloor \left[\alpha A(p)+\left(\beta -\alpha \right)B(p)\right]+\\ \{\begin{array}{ll}\left\{\frac{t}{\beta }\right\}\beta A(p)\hfill & \text{if\hspace{1em}}\left\{\frac{t}{\beta }\right\}\beta \le \alpha \hfill \\ \alpha A(p)+\left(\left\{\frac{t}{\beta }\right\}\beta -\alpha \right)B(p)\hfill & \text{else}\hfill \end{array}\end{array}$

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.

FIGURE 5

Figure 5. Plot of x₁(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)=\varrho sin(\frac{2\pi t}{T}+n\Delta \phi )\mathrm{\text{for}}n=1,\dots ,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:C²(ℝ, ℝ^N) → ℝ as

$D\left[ϵ\right]:={\displaystyle {\int }_0^{T}d}(t,ϵ,\dot{ϵ})dt$

where

$\begin{array}{l}d(t,ϵ,\dot{ϵ}):=d_1(t,ϵ,\dot{ϵ})+w{d}_2(t,\dot{ϵ}),\\ d_1(t,ϵ,\dot{ϵ}):={\displaystyle {\int }_0^{NL}-}\frac{1}{\mu }{f}_{ref}(X,t)v(\chi (X,t))dX,\\ d_2(t,\dot{ϵ}):={\displaystyle \sum _{n=1}^N{\dot{ϵ}}_n^2}(t)\end{array}$

i.e., the energy cost is the time integral over a period of a dissipation rate which is sum of two terms: $d_1(t,ϵ,\stackrel{°}{ϵ})$ is $\frac{1}{\mu }$ times the energy expended to overcome the friction force and $d_2(t,\stackrel{°}{ϵ})$ is the cost of control weighted by a scalar factor w. D[ϵ] is thus $\frac{1}{\mu }$ times the sum of the work due to the friction force plus the L²-norm of the controls suitably weighted to time the input direction.

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

$d_1(t,ϵ,\dot{ϵ})=\dot{ϵ}·\mathbb{D}(ϵ)\dot{ϵ}$

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

$d_2(t,\dot{ϵ})=\dot{ϵ}·{\mathbb{I}}_N\dot{ϵ}$

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

$\begin{array}{cc}D\left[ϵ\right]={\displaystyle {\int }_0^T\dot{ϵ}}·\mathbb{G}(ϵ)\dot{ϵ}dt& (28)\end{array}$

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

Non-Linear Optimal Control Problem

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

$\begin{array}{l}\underset{ϵ\in \mathcal{S}}{\max }\Delta x_0:={\displaystyle {\int }_0^{T}v(ϵ)·\dot{ϵ}}dt\\ \mathcal{S}=\left\{ϵ\in C^2|ϵ(0)=ϵ(T)\wedge D[ϵ]=c\right\}.& (29)\end{array}$

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

$\begin{array}{cc}\frac{d}{dt}\frac{\partial ℱ}{\partial {\dot{ϵ}}_n}(t,ϵ,\dot{ϵ})-\frac{\partial ℱ}{\partial {ϵ}_n}(t,ϵ,\dot{ϵ})=0& (30)\end{array}$

where ${F}(t,ϵ,\stackrel{°}{ϵ}):=\text{v}(ϵ)·\stackrel{°}{ϵ}-\lambda \stackrel{°}{ϵ}·G(ϵ)\stackrel{°}{ϵ}$, λ 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

$\Delta x_0\simeq V[ϵ,\dot{ϵ}]:={\displaystyle {\int }_0^T\mathbb{V}}ϵ·\dot{ϵ}dt.$

and the energy functional by

$\begin{array}{cc}D[ϵ]={\displaystyle {\int }_0^T\dot{ϵ}}·\mathbb{G}(ϵ)\dot{ϵ}dt\simeq {\displaystyle {\int }_0^T\dot{ϵ}}·\mathbb{G}\dot{ϵ}dt& (31)\end{array}$

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

${\left\{\mathbb{G}\right\}}_{ij}=\{\begin{array}{ll}\frac{L^3}{4N}(2i-1)\left(2(N-j)+1\right)\hfill & \text{if }i<j\hfill \\ \frac{L^3}{12N}\left[4N(3i-2)-3{\left(2i-1\right)}^2\right]+w\hfill & \text{if }i=j\hfill \\ \frac{L^3}{4N}(2j-1)\left(2(N-i)+1\right)\hfill & \text{if }i>j\hfill \end{array}$

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

${ϵ}_n^{\star }(t)={\varrho }_a{\varrho }_n\sin \left(\frac{2\pi }{T}t+{\vartheta }_n\right).$

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

$\begin{array}{cc}{\varrho }_{N+1-n}={\varrho }_n\forall n=1,\dots ,N;& (32)\end{array}$

• phase differences between adjacent segments are symmetric about the center, i.e.,

${\vartheta }_{n+1}-{\vartheta }_n={\vartheta }_{N+1-n}-{\vartheta }_{N-n}\forall n=1,\dots ,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).

FIGURE 6

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 e_n for n = 1, …, N can be approximated by a constant average value (cf. Figure 6), i.e.,

${\varrho }_n\simeq \overline{\varrho }\text{\hspace{1em}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 φ₀

$\begin{array}{cc}{\vartheta }_n\simeq n{\vartheta }^{\star }+{\vartheta }_0& (33)\end{array}$

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^{\star }(t)={ϵ}^{\star }(Y_n,t)$

where $Y_n=\frac{X_n+X_{n+1}}{2}$ is the midpoint of the n-th segment and

$\begin{array}{l}{ϵ}^{\star }(X,t)={\varrho }_a\varrho (X)\sin \left(\frac{2\pi }{T}t+\vartheta (X)\right)\\ \simeq {\varrho }_a\overline{\varrho }\sin \left(\frac{2\pi }{T}t+{\vartheta }^{\star }X+{\vartheta }_0\right)\\ ={\varrho }_a\overline{\varrho }\sin \left({\vartheta }^{\star }\left[X-\left(-\frac{1}{{\vartheta }^{\star }}\frac{2\pi }{T}\right)t\right]+{\vartheta }_0\right)\end{array}$

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

$\begin{array}{cc}{ϵ}^{\star }(X,t)\simeq H(X-vt)& (34)\end{array}$

where $H(y):={\varrho }_a\bar{\varrho }sin({\vartheta }^{\star }y+{\vartheta }_0)$ and $v:=-\frac{1}{{\vartheta }^{\star }}\frac{2\pi }{T}$.

FIGURE 7

Figure 7. Plot of piecewise constant optimal strain ${ϵ}_n^{\star }(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).

FIGURE 8

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)=\chi (X,t)=x_0(t)+X+u(X,t).$

The hypothesis of periodicity leads to

$\begin{array}{cc}u(X+NL,t)=u(X,t)\forall X,t.& (35)\end{array}$

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

$\begin{array}{cc}{\displaystyle {\int }_0^{NL}ϵ}(X,t)dX=0\forall t& (36)\end{array}$

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

$\begin{array}{cc}{\displaystyle \sum _{n=1}^N{ϵ}_n}(t)=0\forall t.& (37)\end{array}$

The optimal control problem becomes

$\begin{array}{l}\underset{ϵ\in {\mathcal{S}}^{\prime }}{\max }V[ϵ,\dot{ϵ}]:={\displaystyle {\int }_0^T\dot{ϵ}·\mathbb{V}ϵ}dt\\ {\mathcal{S}}^{\prime }=\{ϵ{{\displaystyle \int }}_0^T\in C^2\left(ℝ,{ℝ}^N\right)|{\displaystyle \sum _{n=1}^N{ϵ}_n}=0\wedge \\ ϵ(0)=ϵ(T)\wedge {\displaystyle {\int }_0^T\dot{ϵ}·\mathbb{G}\dot{ϵ}}dt=c\}& (38)\end{array}$

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

$\begin{array}{cc}e=\left[\begin{array}{c}\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_n\\ ⋮\\ e_N\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{e}_1\end{array}\\ e^{i\frac{2\pi k}{N}}{e}_1\\ ⋮\\ e^{i\frac{2\pi k}{N}(n-1)}{e}_1\\ ⋮\\ e^{i\frac{2\pi k}{N}(N-1)}{e}_1\end{array}\right]& (39)\end{array}$

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

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

- each component can be obtained from the previous one by a rotation of $\frac{2\pi k}{N}$ or, in other words, the phase difference between two consecutive components is constant, i.e., for n = 1, …, N

where ${\vartheta }^{\star }:=\frac{2\pi k}{N}$ and ${\vartheta }_0:=arg(e_1)-\frac{2\pi k}{N}$;

FIGURE 9

Figure 9. 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 u_n through the periodic version of transformation (15), i.e.,

$\begin{array}{cc}ϵ(t)={\mathbb{J}}_{per}u(t)& (40)\end{array}$

where

$\begin{array}{cc}{\mathbb{J}}_{per}:=\frac{1}{L}\left[\begin{array}{cccc}1& & & -1\\ -1& 1& & \\ & \ddots & \ddots & \\ & & -1& 1\end{array}\right].& (41)\end{array}$

In particular, we get

$\begin{array}{l}\underset{ϵ\in {\mathcal{S}}_u^{\star }}{\max }V[u,\dot{u}]:={\displaystyle {\int }_0^T\dot{u}·{\mathbb{V}}_u^{\star }u}dt\\ {\mathcal{S}}_u^{\star }=\{u\in {{\displaystyle \int }}_0^T{C}^3\left(ℝ,{ℝ}^N\right)|u(0)=u(T)\\ \wedge E\left[u,\dot{u}\right]:={\displaystyle {\int }_0^T\dot{u}·{\mathbb{G}}_u^{\star }\dot{u}}dt=c\}& (42)\end{array}$

where

$\begin{array}{l}{\mathbb{V}}_u^{\star }:=J_{per}^T\mathbb{V}{J}_{per}\\ {\mathbb{G}}_u^{\star }:=J_{per}^T\mathbb{G}{J}_{per}\end{array}$

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

$\begin{array}{l}\underset{ϵ\in {\mathcal{S}}_u}{\max }V[u,\dot{u}]:={\displaystyle {\int }_0^T\dot{u}·{\mathbb{V}}_{u}u}dt\\ {\mathcal{S}}_u=\{{{\displaystyle \int }}_0^{T}u\in C^3\left(ℝ,{ℝ}^N\right){{\displaystyle \int }}_0^T|u(0)=u(T)\\ \wedge E\left[u,\dot{u}\right]:={\displaystyle {\int }_0^T\dot{u}·{\mathbb{G}}_u\dot{u}}dt=c\}\end{array}$

where

$\begin{array}{l}{\mathbb{V}}_u:=J^T\mathbb{V}J\\ {\mathbb{G}}_u:=J^T\mathbb{G}J\end{array}$

are “quasi-circulant” matrices indeed

$\begin{array}{l}{\mathbb{V}}_u:={\mathbb{V}}_u^{\star }+{\mathbb{E}}_V\\ {\mathbb{G}}_u:={\mathbb{G}}_u^{\star }+{\mathbb{E}}_G\end{array}$

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

and

where a = 2L³(3−N)+12Nw and b = L³(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[ϵ,\dot{ϵ}]:={\displaystyle {\int }_0^T\left(d_1(t,ϵ,\dot{ϵ})+w{d}_2(t,ϵ,\dot{ϵ})\right)}dt=\overline{c}$

and we let w vary from 0 to 10² 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 $\frac{N}{2\pi }arccos(\frac{1}{2}\frac{6w-L^3}{3w+L^3})$ 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 $\frac{N}{3}$, i.e., one full wave-length every three segments.

FIGURE 10

Figure 10. 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” u_s among harmonic shape functions having the form (in our notation)

$\begin{array}{cc}{ϵ}_n(t)=a\sin \left(\frac{2\pi }{T}t+{\eta }_n\right)\text{\hspace{1em}for }n=1,\dots ,N& (43)\end{array}$

where $a\in (0,\frac{1}{L})$ 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

$\begin{array}{cc}{u}_s=\frac{\Delta x_0}{T}=\frac{1}{T}{\displaystyle {\int }_0^{T}v(ϵ)·\dot{ϵ}}dt,& (44)\end{array}$

the optimization problem reads

$\begin{array}{cc}\underset{\eta \in {\left[0,2\pi \right)}^N}{\max }{u}_s(\eta )=\frac{1}{T}{\displaystyle {\int }_0^{T}v(ϵ)·\dot{ϵ}}dt& (45)\end{array}$

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

$\begin{array}{cc}\underset{\eta \in {\left[0,2\pi \right)}^N}{\max }{\displaystyle {\int }_0^T\dot{ϵ}·\mathbb{V}ϵ}dt.& (46)\end{array}$

Denote the actuation phase differences between adjacent segments by

$p_n:=\Delta \eta (n)={\eta }_{n+1}-{\eta }_n\text{\hspace{1em}for }n=1,\dots ,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

$\begin{array}{cc}{p}_n=p_{N-n}\forall n.& (47)\end{array}$

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

$\widetilde{ϵ}(t):={\left\{{\widetilde{ϵ}}_n(t):=a\sin \left(\frac{2\pi }{T}t+{\tilde{\eta }}_n\right)\right\}}_{n=1,\dots ,N},$

and consider the shape change $\widehat{ϵ}(t)$ associated with

$\widehat{\eta }:=-\mathbb{K}\tilde{\eta }+2\pi$

where

$\mathbb{K}:=\left[\begin{array}{c}\begin{array}{ccccc}0& 0& \cdots & 0& 1\\ 0& 0& \cdots & 1& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \text{0}& \text{1}& \cdots & \text{0}& \text{0}\\ \text{1}& \text{0}& \cdots & \text{0}& \text{0}\end{array}\end{array}\right]\in {ℝ}^{N\times N}.$

Notice that for n = 1, …, N,

${\widehat{ϵ}}_n(t):=a\sin \left(\frac{2\pi }{T}t-{\left(\mathbb{K}\tilde{\eta }\right)}_n\right)={\left\{-\mathbb{K}\widetilde{ϵ}(-t)\right\}}_n$

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

$\begin{array}{l}{\displaystyle {\int }_0^T\stackrel{.}{\widehat{ϵ}}}·\mathbb{V}\widehat{ϵ}dt=-{\displaystyle {\int }_0^T\stackrel{.}{\widetilde{ϵ}}}(-t)·{\mathbb{K}}^T\mathbb{V}\mathbb{K}\widetilde{ϵ}(-t)dt\\ ={\displaystyle {\int }_{-T}^0\stackrel{.}{\widetilde{ϵ}}}·\mathbb{V}\widetilde{ϵ}dt={\displaystyle {\int }_0^T\stackrel{.}{\widetilde{ϵ}}}·\mathbb{V}\widetilde{ϵ}dt.\end{array}$

Thus

$\tilde{\eta }=-\mathbb{K}\tilde{\eta }+2\pi$

which leads to (47).

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

${\displaystyle {\int }_0^T\dot{ϵ}}·\dot{ϵ}dt=\frac{2N}{T}{(a\pi )}^2=:c^{\star }$

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

$\begin{array}{l}\underset{{}_{ϵ\in \mathcal{S}}}{\max }V\left[ϵ,\dot{ϵ}\right]:={\displaystyle {\int }_0^T\dot{ϵ}·\mathbb{V}ϵ}dt\\ \mathcal{S}=\{{{\displaystyle \int }}_0^Tϵ\in C^2\left(ℝ,{ℝ}^N\right)|ϵ(0)=ϵ(T)\wedge {\displaystyle {\int }_0^T\dot{ϵ}·\dot{ϵ}}dt=c^{\star }\}.& (48)\end{array}$

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

$\begin{array}{cc}{ϵ}_n^{\star }(t)=a\sqrt{N}\left|\right|e_n\left|\right|\sin \left(\frac{2\pi }{T}t+\arg \left(e_n\right)+{\vartheta }_a\right)& (49)\end{array}$

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.

FIGURE 11

Figure 11. Average velocities u_s, (44) obtained by the solution $\widetilde{ϵ}$ 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⁻¹⁰, 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.

Funding

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

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.

Acknowledgments

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., and Raff, M. (2002). Molecular Biology of the Cell, 4th Edn. New York, NY: Garland Science.

Google Scholar

Boxerbaum, A., Shaw, K., Chiel, H., and Quinn, R. (2012). Continuous wave peristaltic motion in a robot. Int. J. Robot. Res. 31, 302–318. doi: 10.1177/0278364911432486

CrossRef Full Text | Google Scholar

Boyle, J. H., Berri, S., and Cohen, N. (2012). Gait modulation in c. elegans: an integrated neuromechanical model. Front. Comput. Neurosci. 6:10. doi: 10.3389/fncom.2012.00010

PubMed Abstract | CrossRef Full Text | Google Scholar

Collar, A. (1962). On centrosymmetric and centroskew matrices. Q. J. Mechan. Appl. Math. 15, 265–281. doi: 10.1093/qjmam/15.3.265

CrossRef Full Text | Google Scholar

Daltorio, K., Boxerbaum, A., Horchler, A. S., Shaw, K. M., Chiel, H. J., and Quinn, R. D. (2013). Efficient worm-like locomotion: slip and control of soft-bodied peristaltic robots. Bioinspir. Biomimet. 8:035003. doi: 10.1088/1748-3182/8/3/035003

PubMed Abstract | CrossRef Full Text | Google Scholar

Denny, M. (1980). The role of gastropod pedal mucus in locomotion. Nature 285:160. doi: 10.1038/285160a0

CrossRef Full Text | Google Scholar

DeSimone, A., Guarnieri, F., Noselli, G., and Tatone, A. (2013). Crawlers in viscous environments: linear vs non-linear rheology. Int. J. Non Linear Mech. 56, 142–147. doi: 10.1016/j.ijnonlinmec.2013.02.007

CrossRef Full Text | Google Scholar

DeSimone, A., and 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. doi: 10.1140/epje/i2012-12085-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, C., Hendrix, P., and Arancon, N. (2018). Biology and Ecology of Earthworms. Springer. Available online at: https://www.springer.com/in/book/9780387749426

Google Scholar

Fang, H., Li, S., Wang, K. W., and Xu, J. (2015). Phase coordination and phase–velocity relationship in metameric robot locomotion. Bioinspir. Biomimet. 10:066006. doi: 10.1088/1748-3190/10/6/066006

CrossRef Full Text | Google Scholar

Gardner, C. (1976). The neuronal control of locomotion in the earthworm. Biol. Rev. 51, 25–52. doi: 10.1111/j.1469-185X.1976.tb01119.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Garrey, W., and Moore, A. (1915). Peristalsis and coordination in the earthworm. Am. J. Physiol. Legacy Cont. 39, 139–148. doi: 10.1152/ajplegacy.1915.39.2.139

CrossRef Full Text | Google Scholar

Ge, J., Calderón, A., and Pérez-Arancibia, N. (2017). An earthworm-inspired soft crawling robot controlled by friction. arXiv [preprint] arXiv:1707.04084.

Google Scholar

Gray, J., and Lissmann, H. (1938). Studies in animal locomotion: Vii. locomotory reflexes in the earthworm. J. Exp. Biol. 15, 506–517.

Google Scholar

Grillner, S. (2006). Biological pattern generation: the cellular and computational logic of networks in motion. Neuron 52, 751–766. doi: 10.1016/j.neuron.2006.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ijspeert, A. (2008). Central pattern generators for locomotion control in animals and robots: a review. Neural Netw. 21, 642–653. doi: 10.1016/j.neunet.2008.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Marder, E., and Bucher, D. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, R986–R996. doi: 10.1016/S0960-9822(01)00581-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McInerney, A. (2013). First Steps in Differential Geometry: Riemannian, Contact, Symplectic. Undergraduate Texts in Mathematics. New York, NY: Springer.

Google Scholar

Menciassi, A., Accoto, D., Gorini, S., and Dario, P. (2006). Development of a biomimetic miniature robotic crawler. Auton. Robots 21, 155–163. doi: 10.1007/s10514-006-7846-9

CrossRef Full Text | Google Scholar

Nemitz, M. P., Mihaylov, P., Barraclough, T. W., Ross, D., and Stokes, A. A. (2016). Using voice coils to actuate modular soft robots: wormbot, an example. Soft Robot. 3, 198–204. doi: 10.1089/soro.2016.0009

PubMed Abstract | CrossRef Full Text | Google Scholar

Pehlevan, C., Paoletti, P., and Mahadevan, L. (2016). Integrative neuromechanics of crawling in D. melanogaster larvae. Elife 5:23. doi: 10.7554/eLife.11031

CrossRef Full Text | Google Scholar

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.

PubMed Abstract | Google Scholar

Umedachi, T., Kano, T., Ishiguro, A., and Trimmer, B. (2016). Gait control in a soft robot by sensing interactions with the environment using self-deformation. R. Soc. Open Sci. 3:160766. doi: 10.1098/rsos.160766

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Brunt, B. (2004). The Calculus of Variations. New York, NY: Springer.

Google Scholar

Wang, K., Yan, G., Ma, G., and Ye, D. (2009). An earthworm-like robotic endoscope system for human intestine: design, analysis, and experiment. Anna. Biomed. Eng. 37, 210–221. doi: 10.1007/s10439-008-9597-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, Q., Po, M. D., Hulme, E., Chen, S., Liu, X., Kwok, S. W., et al. (2012). Proprioceptive coupling within motor neurons drives c. elegans forward locomotion. Neuron 76, 750–761. doi: 10.1016/j.neuron.2012.08.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiezel, O., Giraldi, L., DeSimone, A., Or, Y., and 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.

Google Scholar