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.


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

Model Description and Kinematics
Consider a 1D crawler moving along a straight line and assume its reference configuration is the segment [X 1 : = 0, X 2 : = L] , cf.  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., where x 1 (t) : = χ(X 1 , t) and x 2 (t) : = χ(X 2 , t). By definition, s(0, t) ≡ 0 for all t and we assume that in order to guarantee the monotonicity of χ(·, t) at any time t.
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 where we have implicitly defined the strain in terms of which condition (1) reads Finally, notice that, since the material (or Lagrangian) velocity iṡ 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 x 1 (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 where µ > 0 is a friction (or viscosity) coefficient; • or a more general "p-model" 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).
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., Then, by neglecting inertia, the force balance yields 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 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ẋ 1 (t), we obtaiṅ which, in the case of zero external forces, is independent of µ. Notice that, once the initial position x 1 (0) the strain ǫ(X, t) and the external force F e (t) are provided, the whole dynamics x 1 (t) can be determined by integrating (5). Indeed, since and hence the right hand side of (5) is known once ǫ(X, t) is specified.

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.
whence the definition of strain 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., Consequently its time-derivative,ǫ(·, t), is piecewise constant and hence, from (7),ṡ(·, t) is piecewise affine, i.e., for X ∈ [X n−1 , X n ]. 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 where, in view of (9), (1 + ǫ j ) 1−p and, in view of (10), Solving forẋ 0 (t) we obtaiṅ which can be rewritten in the following vectorial forṁ Frontiers in Robotics and AI | www.frontiersin.org Equation (11) fully describes the dynamics once x 0 (0) and ǫ(X, t) are provided. In particular, the displacement after T time units is given by

Relative Displacement
We can rewrite everything in terms of a displacement relative to x 0 (t), defined as in order to describe the displacement in a coordinate system which is "co-moving" with the left end x 0 (t).

Optimal Control Problems
In this section we address the problem of maximizing the net displacement x 0 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 smalldeformation 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 2 function defined from R to R 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ǫ) where A and B are symmetric and positive definite Ndimensional matrices. Overall, the feasible region is

Optimization problem
The general (non-linear) optimization problem is 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,

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., and, integrating by parts, Frontiers in Robotics and AI | www.frontiersin.org In particular, it can be proved that the (skew-symmetric Toeplitz) matrix skw v ǫ (0) = : V depends only on N, L and p. Indeed, Therefore, in the regime of small deformations, problem (18) can be replaced by the following linear problem The corresponding Euler-Lagrange equations lead to the following system of second order linear ODEs In general, a solution to (21) might be difficult to determine due to the complexity of finding a common diagonalization of A and B. However, following the procedure adopted by Wiezel et al. (2018), we can solve this problem when one of the two operators is null, say A ≡ 0 (resp. B ≡ 0), and the other one, B (resp. A), is symmetric, positive definite and such that the eigenspaces associated with the maximum-modulus eigenvalues where α ∈ C \ {0} is a constant such that ||α|| = c 2T and e = (e 1 , e 2 , . . . , e N ) T ∈ C N \ {0} is a suitable constant vector depending only on A and V.
• for A symmetric and positive definite and B = 0, a solution of (20) with ǫ of unitary time frequency must be of the form where α ∈ C \ {0} is a constant such that ||α|| = c 2T and e = (e 1 , e 2 , . . . , e N ) T ∈ C N \ {0} is a suitable constant vector depending only on A and V.
Both (22) and (23) have the form i.e., they are circles in the plane ℜ(e), ℑ(e) , regardless the number of links. Moreover, using the polar representationŝ α = ̺ a e iϑ a and e n = ̺ n e iϑ n ∀n, we get, for n = 1, . . . , N, 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 ̺ n sin 2π T t + ϑ n , thus reducing the number of parameters to 2N + 1. The problem for A = 0 and B = I N is essentially equivalent to the one for A = I N and B = 0 (provided that unitary time frequency of ǫ is prescribed). Indeed, if (25) is a solution to T 2 and vice versa.
In general the two problems, A = 0 with B symmetric positive definite and B = 0 with A 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), Bǫ ·ǫ dt and, analogously,

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 ǫ(X, t) : = ǫ 0 cos 2π L (X − ct) or equivalently, in terms of the stretch, 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, Finally, in view of (5), The Newtonian case is recovered by setting p = 0, i.e.,

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.

Square Contraction Wave
Consider the square contraction wave or equivalently, in terms of the stretch, 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 y mod L).
By integrating the stretch over space, we get Finally, in view of (5), .
Defining α : = L−ξ c and β : = L c , it follows that 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.

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π t T + n ϕ) for n = 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 : C 2 (R, R N ) → R as Frontiers in Robotics and AI | www.frontiersin.org where d(t, ǫ,ǫ) : = d 1 (t, ǫ,ǫ) + w d 2 (t,ǫ), 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, ǫ,ǫ) is 1 µ times the energy expended to overcome the friction force and d 2 (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 L 2 -norm of the controls suitably weighted to time the input direction.
Some calculations [see section 1 in Appendix C (Supplementary Material)] lead to where D(ǫ) ∈ R N×N for any ǫ ∈ (−1, +∞] N , and d 2 (t,ǫ) =ǫ · I Nǫ where I N is the N-dimensional identity matrix. Therefore the energy functional is where G(ǫ) : = D(ǫ) + w I N .

Non-linear Optimal Control Problem
The non-linear optimization problem associated with energy functional (28) is The Euler-Lagrange equations lead to a second order non-linear system of ODEs, i.e., for n = 1, . . . , N,

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 Vǫ ·ǫ dt.
where G : = G(0). Hence, in the small-deformation regime, the problem fits the form (17)-(18) for A = 0 and B = G. Moreover, G is bisymmetric (namely, symmetric about both of its diagonals) and depends only on N, L and w. Indeed Therefore a solution must be of the form (25) ǫ ⋆ n (t) = ̺ a ̺ n sin 2π T t + ϑ n .
The centrosymmetry of G and the skew-centrosymmetry of V 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; Frontiers in Robotics and AI | www.frontiersin.org • 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).
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., 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 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., where Y n = X n +X n+1 2 is the midpoint of the n-th segment and whence the canonical form of traveling wave which describes peristalsis (cf. Figure 7) where H(y) : = ̺ a̺ sin ϑ ⋆ y + ϑ 0 and v : = − 1 ϑ ⋆ 2π T .

The Edge-Effect
The symmetric structure of the optimal gait (in the smalldeformation regime) arises from underlying physical symmetries which clearly stand out in the properties of the matrices G and V. 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 shiftinvariance 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).
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., The hypothesis of periodicity leads to 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 and, in the discrete framework (9), this leads to The optimal control problem becomes for some k ∈ {1, . . . , N − 1} and e 1 ∈ C \ {0}, cf. Figure 9. In particular, -each component of e has modulus ̺ : = ||e 1 ||; -each component can be obtained from the previous one by a rotation of 2π k N or, in other words, the phase difference between two consecutive components is constant, i.e., for n = 1, . . . , N arg(e n ) = (n − 1) 2π k N + arg(e 1 ) = nϑ ⋆ + ϑ 0 where ϑ ⋆ : = 2π k N and ϑ 0 : = arg(e 1 ) − 2π k N ; 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., where In particular, we get where 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 ǫ∈S u V [u,u] where E V and E G are null a part from the last column and the last row, i.e.,

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[ǫ,ǫ] : = T 0 d 1 (t, ǫ,ǫ) + wd 2 (t, ǫ,ǫ) dt =c and we let w vary from 0 to 10 2 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 N 2π arccos 1 2 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 N 3 , i.e., one full wave-length every three segments.
This behavior is qualitatively unaffected by the type of friction model which is adopted (i.e., by the choice of the parameter p).

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) ǫ n (t) = a sin 2π T t + η n for n = 1, . . . , N where a ∈ 0, 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 the optimization problem reads and in the small-deformation regime it can be replaced by Denote the actuation phase differences between adjacent segments by p n : = η(n) = η n+1 − η n for n = 1, . . . , N − 1.
In fact these properties can be rigorously proved under the assumption that problem (45)  where Notice that for n = 1, . . . , N, ǫ n (t) : = a sin 2π T t − (Kη) n = −Kǫ(−t) n and hence, by exploiting the fact that V is skew-centrosymmetric (namely, which leads to (47). Problem (46) constrains the L 2 -norm of the time-derivatives, i.e., for strains having the form (43) we get T 0ǫ ·ǫ dt = 2N T (aπ) 2 = : c ⋆ regardless of η. Therefore we can extend the maximization to the C 2 periodic strains whose time derivative fulfills the same constraint, i.e., (48) Since problem (21) reduces to (48) when A = 0 and B = I N , a solution to (48) must be of the form where e = (e n ) 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.

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.