Abstract
Limbless crawling is a fundamental form of biological locomotion adopted by a wide variety of species, including the amoeba, earthworm and snake. An interesting question from a biomechanics perspective is how limbless crawlers control their flexible bodies in order to realize directional migration. In this paper, we discuss the simple but instructive problem of peristalsislike locomotion driven by elongation–contraction waves that propagate along the body axis, a process frequently observed in slender species such as the earthworm. We show that the basic equation describing this type of locomotion is a linear, onedimensional diffusion equation with a time–spacedependent diffusion coefficient and a source term, both of which express the biological action that drives the locomotion. A perturbation analysis of the equation reveals that adequate control of friction with the substrate on which locomotion occurs is indispensable in order to translate the internal motion (propagation of the elongation–contraction wave) into directional migration. Both the locomotion speed and its direction (relative to the wave propagation) can be changed by the control of friction. The biological relevance of this mechanism is discussed.
1. Introduction
Crawling without legs is one of the most fundamental types of biological locomotion and can be observed in a wide variety of species, including the amoeba, earthworm, snail and leech, as well as in mobile cells [1–6]. The body of a limbless crawler is in contact with the substrate on which it moves over a long distance and continuously interacts with it. This type of locomotion looks qualitatively different from the bipedal or quadrupedal walking of larger animals; an obvious question to ask is how crawlers control their flexible bodies, which have a continuum of degrees of freedom. This question is also related to a current topic in biologically inspired robotics [7], which aspires to create autonomous robots by mimicking animal motions. When crawlers are migrating, forwards–backwards symmetry (in terms of the motion and/or the shape of the body) must be broken in order to select a particular direction of movement. For instance, amoebas extend pseudopods in the direction of movement. For earthworm and leech, elongation–contraction waves that resemble peristaltic motion propagate in a fixed direction along the body axis.
In this paper, we focus on the peristalsislike locomotion typically observed in earth worm. We consider that the category of the peristalsislike locomotion includes the locomotion of snail and of some kind of snakes that can move without undulation of the body (e.g. boas and vipers [3]), although it might be inconsistent with the traditional classification [8,9]. The body structure and size of the above animals are very different from each other: for snail and the abovementioned snakes, the propagation of elongation–contraction wave occurs only around the ventral surface contacting with the substrate. However, we can expect that a common mechanical principle exists in their locomotion, because the internal motion just above the contact interface should be relevant for the locomotion. Hereafter, we inclusively refer the locomotion of the above slender animals as peristaltic locomotion. Also, the locomotion of some kinds of multilegged animals, such as caterpillar and myriapod, has a similarity with peristaltic locomotion (details will be discussed in §6.1).
It has been well recognized that in order to realize peristaltic locomotion, a peristaltic crawler should form anchors (at which the body grips the substrate to increase sliding resistance locally) on some parts of the contact area. For example, the contracted (thickened) body segments of the earthworm project bristles from inside in order to increase the friction coefficient [2,4]. For the abovementioned snakes, slight corrugation of the ventral surface occurs coordinating the peristaltic wave; the parts pressed to the substrate play the role of anchors. A cleverly designed example is the snail locomotion [10–12], which depends on the shear thinning behaviour of the mucus that is secreted on the contact interface. This behaviour reduces the frictional resistance of the parts of the contact interface where a larger shear rate is generated by the peristaltic wave, while the low shear rate parts, where the mucus exhibits relatively high frictional resistance, work as anchors. In the abovementioned multilegged animals, the anchor effect is caused by touching and gripping the substrate with legs.
Interestingly, the direction in which the peristaltic (elongation–contraction) wave propagates depends on the species. For example, when the terrestrial snail (Helix), sea snail (Pomatias) and ear shell (abalone Haliotis) move forward, the peristaltic wave travels in the tailtohead direction, whereas it travels in the headtotail direction when the earthworm, eelworm (Criconemoides), leech [1,2] and the abovementioned snake [3] advance. Pioneering biologists [3,4,9] have diagrammatically analysed the kinematics of peristaltic locomotion, and proposed that the locomotion direction is determined by the phase of peristaltic wave at which the anchors are formed. That is, when the anchor is formed in the contracted parts of the body, the locomotion occurs into the direction opposite to the peristaltic wave propagation. This prediction is consistent with the general tendency of actual peristaltic crawlers. However, as the earlier kinematic analysis relies on diagrammatical and verbal considerations, there remain several ambiguities: for example, the arguments are fairly complicated except a specially simple case of rectangular peristaltic wave; mechanical foundation is lacking.
This paper aims at providing a clear mathematical expression for the relationship between locomotion direction and anchoring phase. For this aim, we propose a simple mechanical model for peristaltic locomotion. The body of a peristaltic crawler is regarded as a onedimensional string located on a substrate exerting a viscous friction force on each part of the system. The sinusoidal peristaltic wave is generated by an active internal tension that mimics muscular action causing peristaltic movement. The anchoring effect is represented by a localized increase of the viscous friction coefficient at the points where the phase of the peristaltic wave has a particular value. With some simplifications (such as neglecting inertial effects), the equation of motion can be given by a linear onedimensional diffusion equation with a source term. The diffusion coefficient is inversely proportional to the viscous friction coefficient, and the source term represents the physiological action that causes the elongation–contraction wave. By a perturbation treatment on the anchor effect, we obtain a mathematical expression that relates the locomotion speed (including the direction) with the anchoring phase and with other physical parameters concerning the peristaltic movement (wavelength and amplitude of peristaltic wave, etc.). Our mechanical model reproduces the prediction of the earlier kinematic theory on the locomotion direction. Based on the obtained result, the energy cost of the peristaltic locomotion is calculated. We also analyse the case where the inertial effect is included. Finally, the applicability and biological implications of our model will be discussed.
2. Importance of friction control for locomotion
2.1. Construction of onedimensional mechanical model to describe the peristaltic locomotion
The aim of this paper is to provide a mechanically and mathematically clear explanation for the relationship between the speed of peristaltic locomotion and the anchoring phase (and other parameters characterizing the internal peristaltic movement). For this aim, we need to construct a mechanical model that is simple enough to allow analytical treatment but captures the essential features of peristaltic locomotion. Here, we explain the outline of our modelling. Figure 1a shows a part of a peristaltic crawler. We imagine a slice region bounded by two virtual cross sections C and C′. The forces exerting on the slice are internal tensions T and T′ on the cross sections, which are produced by endogenous (muscular) action, and a frictional force F_{s} exerted from the substrate. Strictly speaking, endogenous action produces internal stresses distributed over the cross sections; that is not uniform and may have tangential components. For example, bodywall muscles of annelids generate a contraction stress around the periphery of the cross section. For snails and abovementioned snakes, the internal stress should be localized around the venter. In order to consider motion of the slice region into the locomotion direction, however, we may pick up the normal component of the stresses and sum it up over the cross section. The tension T represents the summedup quantity. Representing symbolically the slice by a block, the entire body of the crawler is mechanically equivalent to a train of blocks connected by mechanical elements generating and transmitting tension T (figure 1b).
We choose a type of T that produces a sinusoidal peristaltic wave. The choice of the friction law to determine F_{s} is a thorny problem. If we intend to describe locomotion of a particular peristaltic crawler in detail, we need to find out a complicated friction law which depends on the texture of the ventral surface and on the structure of the contact interface (e.g. the intermediate material is soil particles, viscous mucus or moist grass, etc.). As mentioned above, this paper focuses on the relationship between the phase of peristaltic wave at which anchors are formed and the direction of locomotion. In our standpoint, the choice of the friction law for unanchored parts is not crucial; it is sufficient that the frictional resistance is small enough compared with the anchored parts. We employ the linear viscous friction law, i.e. F_{s} = −ζv_{s} (v_{s} is the sliding velocity of the slice region; figure 1a). The crucial problem is how to treat the anchoring effect. For real peristaltic crawlers, the anchoring is realized by some gripping devices of which detailed mechanics depends on species; what they have in common is that anchored parts of the body never or hardly slide along the substrate even though the imbalance of the internal tension T − T′ is large. In our viscous friction model, the anchoring effect can be represented by intensively increasing the coefficient ζ of particular parts (slices or blocks) of the system. The simplicity of the treatment is the major advantage of our modelling.
2.2. A case study using an elementary twoblock model
Before going into a detailed treatment on the multiblock model shown in figure 1b, we first consider a very simple twoblock model (figure 1c) in order to demonstrate some essential features of crawling locomotion under noninertial conditions. Analysis of this simple model reveals that (i) adequate control of friction with substrate is indispensable in order for active internal deformation to cause locomotion, and that (ii) if inertial effects are negligible, locomotion occurs under the condition that the total external force exerted on the system is always zero.
We consider a system comprising a pair of rectangular blocks (with a mass m) connected by an active mechanical element (figure 1c). The active element (hereafter, referred to as the ‘active spring’) consists of two linear springs with a spring constant 2K and a hinge whose opening angle can actively be varied. If the opening angle is fixed, the active spring behaves like a normal linear spring with a spring constant K, but we are considering the case where the natural length is temporally modified by the active change of the opening angle of the hinge. We assume that when the opening angle is varied periodically, the natural length of the active spring varies in sinusoidal fashion with respect to time and is given by ℓ_{0} + a sin ωt, where ℓ_{0} is the average length and a is the amplitude of oscillation. The system is placed on a flat and wet substrate. The term ‘wet’ implies that the substrate exerts no dry friction force on the blocks, but instead exerts a dynamical (linear viscous) friction force that is proportional to the sliding velocity of the blocks.
The question that we address here is whether the system can undergo directional motion, or what conditions are required to realize locomotion. The equations of motion are 2.1where X_{n} and ζ_{n} are the position and the viscous friction coefficient of block n, respectively (n = 1 denotes the lefthand block and n = 2 denotes the righthand block). Below, we make the overdamped approximation^{1} for the above equations of motion (i.e. setting ), and treat the simplified equations of motion, , with initial conditions of X_{1} (0) = 0 and X_{2} (0) = ℓ_{0}.
If ζ_{1} = ζ_{2}, the velocity of the centre of the blocks, , is always zero, as proven by summing both sides of the above equations of motion. This fact is obvious from the symmetry of the system.
We now address the case of ζ_{1} ≠ ζ_{2}. Introducing the relative displacement y ≡ X_{2} − X_{1} − ℓ_{0}, we obtain , where τ^{−1} = K(ζ_{1}^{−1} + ζ_{2}^{−1}). The solution is y = p cos(ωt + ϕ) + q exp(−t/τ), where , ϕ = tan^{−1}(1/τω) and q = aτω/(1 + τ^{2} ω^{2}). Thus, is almost sinusoidal (except for the exponentially decaying component), and hence the centre position oscillates forwards and backwards rather than undergoing locomotion.
Intuitively, because the interaction force exerted on each block (∝ K(y − a sin ωt)) is almost sinusoidal, the blocks both undergo forwards–backwards oscillation. The asymmetry of the system (ζ_{1} ≠ ζ_{2}) is reflected by the difference in oscillation amplitudes (the ratio of which is equal to ζ_{1}/ζ_{2}), but it never causes continual motion in one direction.
One may ask what modifications are required to realize directional motion. The simplest (and almost unique) answer is to switch the amplitudes of ζ_{1} and ζ_{2} in response to the sign of the spring force (assuming, for example, that the asperity of the substrate is variable). During the periods in which the active spring generates a repulsive force on the blocks, we tune the friction coefficients such that ζ_{1} ≫ ζ_{2}; block 2 is then pushed forwards. During the periods of attractive force, we make ζ_{1} ≪ ζ_{2}; block 1 is then pulled towards block 2. Such a mode of locomotion is essentially that exhibited by the inchworm.
We next focus on the total external force exerted on the system. Because the system consists of the two blocks and the active spring, the frictional forces from the substrate to the blocks, ζ_{1}Ẋ_{1} and ζ_{2}Ẋ_{2}, should be regarded as external forces. The sum of these external forces is always zero under the overdamped condition. The directional motion of the two block system is caused by the difference in mobility of the blocks (ζ_{1}≠ζ_{2}).
3. Multiblock model for a slender limbless crawler
3.1. Derivation of diffusion equation as equation of locomotion
In order to model the locomotion caused by the propagation of an elongation–contraction wave, we extend the elementary model described in §2 by forming a onedimensional train of blocks and active springs (figure 1b). As in the case of the twoblock model, we assume that the substrate is wet (i.e. it exerts a viscous friction force only) and that all inertial effects are negligible (the inertial effect will be considered in §6.4). The active springs have a common natural length ℓ_{0} in the reference state of the system and a common spring constant K. In order to mimic the biological action that causes the elongation–contraction waves propagating along the body (referred to as ‘endogenous action’), we modulate the natural length of each spring. The equations of motion for our model are 3.1and 3.2where X_{n}(t) and ζ_{n}(t) are the position and coefficient of viscous friction of the nth block, respectively, and ℓ_{n+1/2}(t) is the modified natural length of the active spring between the nth and (n + 1)th blocks. a(<ℓ_{0}) and N are the amplitude and the periodicity (about the index n) of the modification. The dependence of ζ on n and t allows spatiotemporal control of the viscous friction. Introducing the displacement of nth block from its ‘reference position’ nℓ_{0}, U_{n}(t) = X_{n}(t) − nℓ_{0}, the above equations of motion can be rewritten as 3.3
To move on to a continuum description, we introduce a reference coordinate s = nℓ_{0} and take the limit of ℓ_{0} → 0 by making the following quantities remain finite: η(s,t) = ζ_{n}(t)/ℓ_{0}, κ = ℓ_{0}K, e_{0} = a/ℓ_{0} and k = 2π /(Nℓ_{0}). Representing the displacement field in the continuum description as u(s,t) and replacing the difference with respect to n by ℓ_{0}∂_{s}, equation (3.3) yields 3.4and 3.5where the subscript of u denotes the partial derivative with respect to that quantity; and e_{0} (dimensionless), ω and k are the amplitude, frequency and wavenumber of e. There is the restriction that e_{0} < 1, because e_{0} corresponds to the maximum of a/ℓ_{0} and e_{0} ≥ 1 would imply overlap and mutual permeation of the blocks. The physical meaning of equation (3.4) can be clear by writing the righthand side as ∂_{s}κ (u_{s}−e). That is, the difference between the actual strain u_{s} and the spontaneous strain e(s,t) generates internal tension T = κ(u_{s} − e), and spatial inhomogeneity of the tension T_{s} balances with the viscous friction force ηu_{t}. We note that when η(s,t) is a given function, the governing equation of our system, equations (3.4) and (3.5), is a linear diffusion equation with a nonuniform diffusion coefficient D(s,t) = κ/η(s,t) and a travellingwavetype source term arising from the endogenous action.
Below, we assume that the system size is infinite and u(s,t) has the spatial periodicity with the wavelength λ = 2π/k.
First, we consider the case where the friction coefficient is given by a constant η_{0}; from equations (3.4) and (3.5), 3.6where 3.7
The travellingwave solution for this diffusion equation is 3.8 3.9and 3.10where 3.11is the frequency normalized by the diffusion time k^{2}D_{0}. Equation (3.8) shows that the displacement of a material point s,u(s,t), merely oscillates around zero; there is no directional motion of the system as a whole.^{2} In other words, adequate spatial modulation of η is indispensable in order for the propagation of the elongation–contraction wave to cause directional motion.
4. Perturbative treatment of the effect of anchoring on locomotion
4.1. Spatiotemporal control of friction in response to phase of elongation–contraction wave
As mentioned in §3.1, directional motion of the system as a whole never occurs when η is spatially uniform. We will now examine the effect of introducing anchor points on the motion of the system. In order to include the effect of anchoring, we locally increase the viscous friction coefficient; this is modelled using an array of rectangular pulses in the η–s relation, as shown in figure 2a. The spacing of adjacent pulses is equal to the wavelength λ = 2π/k of e, and each pulse moves at the phase velocity given by v = ω/k.
We consider the case in which the amplitudes of the pulses are much smaller than the background level of η_{0}, and thus treat them as a perturbation. By setting η(s,t) = η_{0}(1 + ɛp(s,t)) in equation (3.4), equation (3.6) is modified to 4.1where p(s,t) is a function representing the array of pulses and ɛ is a small parameter (both are dimensionless). Expanding the solution of equation (4.1) in the form u = u^{(0)} + ɛu^{(1)} + … , we obtain a series of diffusion equations with source terms: 4.2and 4.3
The solution of the lowest order equation has already been obtained (see equations (3.8)–(3.11)), which does not give rise to directional motion of the system. The question to be answered here is whether u^{(1)} leads to directional motion. The answer can be obtained by calculating the spatial average of u_{t}^{(1)}(s,t) over the wavelength λ. Because we can expect u^{(1)} has a spatial period of λ (this can be shown by the rigorous treatment of equation (4.3) given in appendix A), averaging both sides of equation (4.3) yields 4.4where we have introduced the following abbreviation: 4.5
We note that 〈u_{t}^{(1)}〉_{λ} is independent of t and of arbitrary shifts of the integral interval. This is because both u^{(0)}(s,t) and p(s,t) have the spatial period λ and move at a constant phase velocity v.
In order to simplify the calculation, we employ an array of δfunctions for p(s,t): 4.6 4.7where B is the shift of a δfunction from a peak in u^{(0)} (figure 2b) and β = kB is the phase shift corresponding to B. The prefactor of 2π has been inserted to satisfy the normalization condition of 〈p(s,t)〉_{λ} = 1. By substituting equations (3.8) and (4.7) into the righthand side of equation (4.4) (it is convenient to set t = α/ω before the integration is carried out), we obtain 4.8
Thus, of the order of O(ɛ^{1}), the average velocity of the system V = ɛ〈u_{t}^{(1)}〉_{λ} is given by 4.9
For obtaining the second equality above, we have used equation (3.10) and ω /k = v (the phase velocity). Because e_{0}, cos α and sin β are all less than unity (and ɛ is a small parameter), the righthand expression for V in equation (4.9) indicates that the locomotion speed cannot exceed the phase velocity v. Note that the limitation of V < v does not come from the perturbative nature of our treatment. The same limitation exists for a nonperturbative variant of the present model introduced in the discussion part (see equation (6.5)).
Owing to the perturbative nature of the above calculation, equation (4.9) is an approximative relationship which is correct for infinitesimally small ɛ. By contrast, in the limit where the rigidity of the system is infinite (κ → ∞) and the endogenous action generating the elongation–contraction wave is strong enough, a relationship that is essentially equivalent to equation (4.9) holds as an exact result; details are provided in §6.2.
4.2. Relationship between locomotion direction and anchoring phase
Equation (4.9) states that if sin β < 0, then V > 0, i.e. the system moves in the direction of propagation of the elongation–contraction wave. The system moves in the direction opposite to the wave propagation (V < 0) when sin β > 0. In order to illustrate the physical meaning of this condition, we show in figure 2b the shapes of p and u^{(0)} (the principal part of the deformation field) over the interval 0 ≤ s ≤ λ where we set t = α/ω such that a peak in u^{(0)} is located at s = 0. Here, B = β/k represents the shift of the impulse position from the position of the peak in u^{(0)}.
The condition that sin β < 0 for V > 0 indicates that the anchor position B exists in the interval λ /2 < s < λ in figure 2b, where u_{t}^{(0)}(B,α/ω) = ωAsin β < 0. Namely, in order for locomotion to be realized along the direction of wave propagation, the anchors should be located in parts of the system where the local velocity (u_{t}^{(0)}) is negative. Likewise, locomotion in the direction opposite to the wave propagation is obtained when the anchors are located in parts with positive local velocity, corresponding to the interval 0 < s < λ/2 in figure 2b.
The reason for the above criteria can be intuitively understood as follows. In the nonanchored (unperturbed) state, the local velocity field u_{t} is selfbalancing: some parts of the system have positive u_{t} and others have negative u_{t}. The average velocity satisfies V = 〈u_{t}〉_{λ} = 0 for any time t. The role of the anchors is to break this balance. Consequently, the entire system moves in the direction of the dominant velocity component. The maximum value of V is accomplished at B = λ/4 (β = π/2) for V < 0 and 3λ /4 (β = 3π /2) for V > 0.
The condition of equation (4.9) can also be interpreted in terms of the local deformation (elongation or contraction) of the system. For the interval λ/2 < s < λ (0 < s < λ/2), the system is elongated (contracted) because the elongation strain u_{s} is positive (negative). We then obtain the following criterion: for crawlers moving with (against) the direction of wave propagation, anchors should be placed at elongated (contracted) parts of the body. Figure 3a shows this criterion schematically. The upper (lower) organism moves with (against) the wave propagation direction by anchoring at elongated (contracted) parts of its body. This is consistent with the result of the kinematic consideration by pioneering biologists [3,4,9].
According to the criterion shown in figures 2b and 3a, the direction of locomotion is reversed if the peristaltic wave changes direction while the anchoring phase is retained. Some real peristaltic (and similar kinds of) crawlers might exhibit this behaviour when reversing their direction of movement. Indeed, we have observed in a preliminary experiment that the millipede reverses the direction of the density wave of its feet in order to go backwards. It will be interesting to test, in a wide range of crawling organisms, whether and how they reverse their direction of movement by tread force measurement [11,13,14].
5. Characterization of locomotion
5.1. Tradeoff between speed and mileage
Thus far, we have focused on the relationship between the mobility of the system and the anchoring position. In this section, we evaluate quality of the locomotion that is generated by optimal anchoring (i.e. when sin β = 1). We are interested in both the average velocity V and a mileage quantity M (defined below). The latter is related to how far the system can move per unit of energy input.
By setting sin β = 1, equation (4.9) yields 5.1
We note equations (3.9) and (3.11); and thus, and . In order to obtain a good measure of the mileage, we must look at the energy budget of the system. Integration of the elastic energy density (κ/2)(u_{s} − e)^{2} over a distance λ yields a constant, because both u and e have the same wavelength λ and the same phase velocity v: By partial integration, we obtain 5.2
This relationship represents the balance between the viscous dissipation rate (the lefthand side) and the energy input rate (the righthand side) resulting from the active change of e (note that κ(u_{s} − e) is the internal tension). Thus, the integrand on the righthand side, denoted by , represents the power density (J (s m)^{−1}). The spatial average of is calculated as 5.3
In the lefthand (approximate) equality above, we set u_{s} ≈ u_{s}^{(0)} ignoring errors of the order of O(ɛ^{1}). We introduce a ‘mileage’ quantity M (m (J m^{−1})^{−1}) as follows: 5.4
Here, M represents the distance moved for one unit of energy density,^{3} thus large M implies low energetic cost for movement.
From equations (5.1) and (5.4), we obtain the result that 5.5
This indicates an inversely proportional relationship between V and M: no matter how the variables e_{0}, k and ω are adjusted, the time and energy efficiencies cannot be improved simultaneously.
5.2. Matching between diffusion dynamics and endogenous action
We now investigate the behaviour of V and M as functions of e_{0}, k and ω, by fixing κ and η_{0} (and thus D_{0}). Note that e_{0}, ω and k characterize the endogenous action of e(s,t) (see equation (3.5)), whereas κ and η_{0} can be regarded as material and environmental parameters.
Figure 4 shows plots of V and M versus e_{0}, ω and k (see equations (5.1) and (5.4)). In each plot, all other variables are set to unity. The following three results are apparent.

— Dependence on e_{0}: we see that V ∝ e_{0} and M ∝ 1/e_{0} (figure 4a). The former relationship reflects the fact that we study a linear system that transforms an input function e(s,t) into an output function u(s,t), thus the amplitude of u is proportional to that of e. The e_{0} dependence of V is inherited from u. The latter relationship arises because (while V ∝ e_{0}). The divergence of M for e_{0} → 0 implies that if the system is allowed to spend an infinitely long time (V → 0 as e_{0} → 0), it can move with an infinitely small energy cost. This nondissipative behaviour for e_{0} → 0 can be regarded as a quasistatic limit.

— Dependence on ω: we see that , namely, that V ∝ ω for small ω, but V saturates at an upper limit for large values of ω (figure 4b). The crossover between these two types of behaviour occurs at approximately , i.e. when ω = k^{2}D_{0}. This indicates that the frequency ω should be matched to (the inverse of) the diffusion time k^{2}D_{0} in order to achieve high time efficiency, but too high a frequency ω is of no use to the organism. As in the case of the e_{0} dependence of M, the divergence of M for ω → 0 corresponds to a quasistatic limit.

— Dependence on . As seen in figure 4c, the relationship between V and k reaches a maximum when (D_{0}k^{2})/ω = 1 is satisfied, corresponding to matching between the wavelength (∝ 1/k) and the diffusion length ().
Figure 4d shows a threedimensional plot of V as a function of ω and k (all other parameters are set to unity). The ‘ridge’ corresponds to the condition of ω = k^{2} (the matching condition for D_{0} = 1). For comparison, we also provide a threedimensional plot of V versus ω and k for the case where the inertial effect is involved (figure 4e; see §6.3).
We note that owing to the inversely proportional relationship between V and M, the most timesaving mode of locomotion is also the most energyconsuming mode (the minimum of M in figure 4c). That is to say, in this case it is true that time is money.
6. Discussion
6.1. Criterion to determine direction of peristaltic locomotion
The essential point of our theoretical analysis is summarized in figure 2b, from which we can extract the following twostep criterion for the direction of peristaltic locomotion. In the first step, we consider a virtual anchorless state. The peristaltic movement gives rise to a periodic velocity field along the system, which corresponds to u_{t}^{(0)}(s,t) in the diffusion equation model. In the second step, we introduce anchors. If the anchors are placed at the points that had a negative velocity in the anchorless state (the term ‘negative’ represents being opposite to the wave propagation), the locomotion occurs into the positive direction, because the anchors cut down the negative velocity component. Another picture for the anchoring effect is that the anchored points pull/push neighbouring parts into the positive direction. This picture is clear and reasonable from a mechanical viewpoint; we can expect that the twostep criterion qualitatively applies to wider cases, regardless of details of friction law. To demonstrate the generality of the criterion, we consider the locomotion of myriapods [3] below.
Figure 5 schematically represents side views of a myriapod. Figure 5a shows a virtual anchorless state where the body of the myriapod is lifted up from the ground. We assume that the swing angle ϕ, defined in the right illustration in figure 5a, changes in the form of a sinusoidal wave travelling rightward: ϕ(n,t) = ϕ_{0} sin(ωt − qn), where n is the number assigned to legs and q is the wavenumber with respect to n. We introduce the lateral displacement of the foot u(n,t) = L_{0} sin ϕ, as shown in the right illustration (L_{0} is the length of the legs in the side view). u(n,t) corresponds to u^{(0)}(s,t) in the viscous friction model. For legs P and Q, u(n,t) = 0. By adequately choosing the origin of t, we can set (ωt − qn) = 0 for P and (ωt − qn) = −π for Q; thus, the lateral velocity of the foot v_{foot} = u_{t}(n,t) = ωL_{0}ϕ_{0} cos ϕ sin(ωt − qn) is positive for P and negative for Q. Then, we consider the anchor. If the foot of Q contacts the ground with nonslip condition while the foot of P is lifted up from the ground by changing the joint angles of the leg (see the right illustration representing a crosssectional view), the main body is driven into the positive direction, opposite to v_{foot} of the leg Q in the anchorless state. Note that the feet of P and of Q exist in regions where the feet density is high and low, respectively. The argument here is essentially identical to what figure 2b explains.
Here, we insert a brief comment on the anisotropic friction. The importance of anisotropic friction is emphasized in some literatures [4,15–18]. For example, bristles of earthworm tilt to increase frictional resistance to backward sliding. Indeed, when the friction is anisotropic, internal peristaltic movement is automatically converted into directional locomotion without active anchor control. This is because of the following reasons: in the virtual anchorless state (e.g. the state where the system is lifted up from the substrate), the internal peristaltic movement makes some parts move into the direction of the stronger friction. When the system returns onto the substrate, the parts automatically serve as anchors owing to the stronger friction. According to our criterion, the system moves into the direction of weak friction. A shortcoming of this manner of locomotion is that the direction of locomotion is fixed from the beginning, and thus the crawler cannot undergo backward movement even though it reverses the direction of the peristaltic wave. At the best of author's knowledge, it is not clear where there exists any real peristaltic crawler that undergoes directional locomotion using only anisotropic friction (without the active control of the friction). For the case of earthworm, the anchor position is actively changed by cycles of poking and retracting the bristles. At least, if a peristaltic crawler undergoes an active control of the anchor position coordinated with peristaltic wave propagation, the friction does not have to be anisotropic.
6.2. Consideration on rigid machinery limit
Equation (4.9), the criterion that determines the direction of locomotion, is an approximative relationship of the order of O(ɛ). For the relationship to apply with high accuracy, the anchoring intensity (or perturbation parameter) ɛ must be so small that the anchor does not disturb the internal deformation. In contrast, if the springs in figure 1b are replaced by rigid rods and the action that changes the opening angle of the hinges is powerful enough, the internal deformation of the system (i.e. the relative displacement between adjacent blocks) is not affected by the frictional force on the blocks. In this ‘rigid machinery’ limit, we expect that a relationship essentially equivalent to equation (4.9) will hold as an exact condition, regardless of the intensity of the anchoring. Below, we consider the rigid machinery limit, and show that this expectation is correct.
In the continuous description, the rigid machinery limit is characterized by κ → ∞ and u_{s} = e(s,t). This implies that the strain field u_{s} can follow e with no phase delay. Instead of equation (3.6), we have 6.1where the tension T(s,t) is an undetermined function. We employ equation (3.5) and the following expression for η(s,t): 6.2and 6.3
We note that q(s,t) does not include the parameter α indicating a phase delay between e(s,t) and u_{s}(s,t), and that the intensity of anchoring, μ > 0, is arbitrary. From u_{s} = e = e_{0} sin(ω t − ks), 6.4where U(t) represents the component of uniform translation. Because we can assume that T(s,t) has a period of λ, the spatial average of both sides of equation (6.1) gives 〈ηu_{t}〉_{λ} = 0. By substituting equations (6.2)–(6.4) into equation 〈ηu_{t}〉_{λ} = 0 and rearranging the terms, we obtain U′(t) = (e_{0}ω/k)〈η(s,t) sin(ωt − ks)〉_{λ}/〈η(s,t)〉_{λ}. Because η(s,t) = η(s − vt, 0) and η(s, 0) has the periodicity λ, we can set t = 0 in the above spatial averages, which indicates that the translational velocity U′(t) is independent of t (hence can be denoted by V). By performing the spatial integration of the average, we obtain 6.5
Here, V increases with increasing anchoring intensity μ, and saturates at an upper limit of e_{0}v sin β as μ → ∞. Because e_{0} < 1 (see the discussion below equation (3.5)), the locomotion speed V cannot be faster than the phase velocity v.
6.3. Applicability of the model to different friction modulations
We have investigated the effect of localized anchoring on directional locomotion by formulating a linear perturbation model expressed in equations (4.1)–(4.7). Owing to the linearity of our model, the results obtained in §4 can be applied to different types of modulation of the friction coefficient. If the sign of the friction pulse is flipped such that p(s,t) → −p(s,t), the sign of V (and thus the criterion for the locomotion direction) is reversed. This follows from the fact that 〈u^{(1)}〉_{λ} depends linearly on p(s,t) (see equation (4.4)). The locomotion expressed by the negative pluses can be described in the following manner. Almost all parts of the body are anchored, but the anchor is locally weakened at the points where the negative pulses are located, which allows faster sliding of the points. The directional motion of the system as a whole is determined by the localized sliding of these released points.
Thus far, we have assumed that only one pulse exists in an interval of the length λ in sspace. Our calculations can be extended to the case where the interval includes m pulses; the friction coefficient is then given by η(s,t) = η_{0}(1 + ɛP(s,t)), where (here we have explicitly represented the βdependence of p). The factor f_{j} (which can be negative) is the intensity of the pulse with phase shift β_{j}(0 ≤ β_{j} < 2π). Because of the linearity, the average velocity in this multipulse case is , where V(β ) = −ɛωA sin β is the average velocity in the single pulse case (equation (4.9)).
An analogous calculation is possible for the case in which the perturbation P is given by a continuous function with spatial period λ and phase velocity v = ω/k, that is, P = P(s,t;α) = f(ks − (ωt − α)), where f(·) is a 2πperiodic function. Noticing that f(ks − (ωt − α)) = (2π)^{−1} ∫_{0}^{2π} dβ f(β) p(s,t;β), we can calculate the average velocity V_{f} by . For example, if f(β) = sin(β − γ), we obtain V_{f} ∝ (−1)cosγ: the direction of locomotion is determined by the phase shift γ that indicates the extent of overlap between the two sine functions f(β) and sin β. We note that V_{f} > 0 if γ satisfies the condition π/2 <γ < 3π/2 (i.e. if a peak in f(ks) = sin(ks − γ) exists in the interval λ /2 < s < λ). This condition corresponds to that of the single pulse case in figure 2b, that is, V > 0 if the δpulse exists in the interval λ/2 < s < λ.
The idea of perturbative treatment for friction modulation can apply to different situations. For example, it is interesting to ask how the randomness of the environment (i.e. the condition of the substrate) affects locomotion. This problem can be treated by including a random position function R(X) = R(s + u) in the perturbation function, such that P(s,t) = f(ks − (ωt − α)) + R(s + u). Here we impose, for example, the white noise condition of 〈R(y)〉_{s} = 0 and 〈R(y)R(y + z)〉_{s} ∝ δ(z) (the subscript ‘s’ represents the statistical (or longdistance spatial) average). The statistical average of 〈u_{t}^{(1)}〉_{s} determines the mean velocity of locomotion.^{4}
6.4. Consideration of inertial and dry friction effects
A major simplification of our model is that we have neglected inertial effects and dry friction. We now briefly discuss these factors. Inertial effects can be included in our model by including the inertial term ρu_{tt} (ρ is the line density) in the lefthand side of equation (3.4). When η(s,t) includes a weak modulation, the governing equation for the system is 6.6
The system has two characteristic timescales of (the inverse of) the diffusion time ω_{d} and the eigenfrequency ω_{c} defined by 6.7where is the velocity of the elastic wave. For convenience of the description below, we introduce a dimensionless parameter 6.8
The perturbation expansion of equation (6.6) yields 6.9and 6.10
Since u_{s}^{(1)} is expected to have the same spatial period λ as e(s,t) and u^{(0)}, the average velocity V can be calculated in the manner completely parallel to the case without the inertial term. The travellingwave solution of u^{(0)} is 6.11 6.12and 6.13
Averaging both sides of equation (6.10) yields 〈u_{t}^{(1)}〉_{λ} =−〈p(s,t)u_{t}^{(0)}(s,t)〉_{λ}, because we have (i.e. the λperiodicity of u_{t}^{(1)} provides a constant value of 〈u_{t}^{(1)}〉_{λ}). When p is given by an array of δfunction, we obtain 〈u_{t}^{(1)}〉_{λ} = −ωAsin β and V = −ɛ ωAsin β. For optimal anchoring ( sin β  = 1), V is given by 6.14
If we regard V as a function of ω fixing other parameters, it reaches a maximum of 6.15for 6.16
We note that equation (6.16) is the condition for occurrence of the resonance between the elastic oscillation and the endogenous action. Equation (6.15) shows that the maximum velocity is proportional to the diffusion coefficient D_{0}. Figure 4e shows a threedimensional plot of V as a function of ω and k (all other parameters including r are set to unity). In the present case, including the inertial effect, both a cross section of V (ω, k) at a constant k and that at a constant ω have a maximum.
In contrast, dry friction is a strongly nonlinear effect. The addition of a new term describing dry friction violates the linearity of our equations of motion and makes them much harder to deal with. In other words, the behaviour of the system becomes more complex, especially, if inertial effects are also included (e.g. a new oscillatory mode caused by stick–slip instability should be involved). Although detailed analysis of the nonlinear equations of motion is beyond the scope of this paper, it is at least possible and informative to make a comparison with the simple kinematic model discussed by biologists [5,9], in which dry friction dominates the locomotion dynamics.
In this model, the body of a peristaltic crawler is represented by a onedimensional array of blocks (figure 3b). Each block (corresponding to a body segment) adopts either an elongated or contracted state; contracted blocks become thicker and are prevented from sliding along the ground by dry friction, whereas elongated (thin) blocks are lifted up from the ground and are able to move freely in the lateral direction. The body as a whole is assumed to adopt a periodic domain structure consisting of alternating sequences of N_{c} consecutive contracted blocks and N_{e} consecutive elongated blocks. When the elongation–contraction wave propagates towards the right, the system as a whole moves in the opposite direction, as seen in figure 3b. Although the basic picture of the kinematic model is fairly different from ours, the same prediction for the locomotion direction is obtained: anchoring contracted parts of the body gives rise to locomotion opposite to the direction of wave propagation. Therefore, we expect that although our model is based on a particular (simple) friction law, the qualitative aspects of the predictions that it yields are generic (see §6.1).
7. Concluding remarks
We close this paper with the following concluding remarks. In the first paragraph of §1, we stated that peristaltic crawlers break forwards–backwards symmetry. Here, we emphasize that much attention should be paid to the character of the symmetry breaking. In the twoblock model with no anchor control and in the continuum model with uniform η, backwards–forwards symmetry is clearly broken (i.e. for an arrangement of two blocks with different viscous friction coefficients and for the propagation of an elongation–contraction wave). However, no directional motion occurs. On the other hand, anisotropic friction seems somewhat excessive; it inhibits backward motion (when the same head direction is maintained). Most real crawling organisms seem to adopt a more moderate process involving the active spatiotemporal control of friction. In this paper, we have considered an essential aspect of the friction control, based on a simple mechanical model describing onedimensional locomotion driven by peristaltic wave propagation. From analytical treatments of the proposed model, we extracted a twostep criterion for locomotion direction. In the first step, we consider a virtual anchorless (or frictionless) state, and look at the distribution of local velocity produced by internal motion. In the second step, we make anchors at parts of which local velocities (in the anchorless state) point a particular direction; the formation of the anchors results in locomotion into the opposite direction. The basic idea of this criterion would apply to various cases other than the onedimensional peristaltic locomotion, for example, the locomotion of ear shell for which the twodimensionality of contact area cannot be ignored [1,3]. We expect the insight obtained in this paper will be helpful to understand a wide class of limbless locomotion.
Acknowledgements
Y.T. thanks M. Kuhara for showing his numerical results at the earlier stage of this study. This research was funded by MEXT KAKENHI 20300105. Also, it is partially supported by Human Frontier Science Programme RGP51/2007 and Japan–Sweden Strategic International Cooperative Programme.
Appendix A
As a result of the perturbative treatment of equation (4.1), we encounter a series of diffusion equation with a diffusion constant D_{0} and source term. The (formal) solution of f_{t} = D_{0} f_{ss} + g(s,t) can be obtained by the following formula: A 1where A 2is the Gaussian Green's function. In equation (4.3), f = u^{(1)} and g = −pu_{t}^{(0)}. Substituting equations (3.8) and (4.6) into this formula, we obtain A 3where Δu^{(1)} = u^{(1)}(s,t) − u^{(1)} (s,0), C = −λωA sin(kβ) and B_{n} = B + nλ − α/k. We note that owing to the infinite sum with respect to n, u^{(1)}(s,t) has the spatial period of λ if u^{(1)}(s,0) has the same special period.
Since for any f(s), the spatial average of Δu^{(1)} is given by A 4
The integration about s yields unity because the integrand is a normalized Gaussian function. Thus, we reproduce
Footnotes

↵1 For the overdamped approximation to be acceptable, a condition should be imposed on the amplitude of the inertial and frictional forces: the order of these forces is estimated as maω^{2} and ζaω (ζ = min{ζ_{1},ζ_{2}}), respectively. So, it is required that maω^{2} ≪ ζaω, that is, 1 ≪ ζ/(mω). When the blocks are cubes with a linear dimension of L and a density ρ, we have ζ = η_{vis} L^{2}/d and m = ρL^{3} (η_{vis} is the viscosity of the liquid on the substrate and d is the gap distance between the substrate and the bottom of the blocks). Thus, the required condition is 1 ≪ η_{vis}/(ρLdω). For example, taking a set of parameter values such as η_{vis} = 10^{−3} Kg m s^{−1} (close to the viscosity of water), ρ = 10^{3} Kg m^{−3} (the density of water), L = 1cm = 10^{−2} m, d = 1 µm = 10^{−6} m and ω = 1 s^{−1}, the condition is satisfied (η_{vis} /(ρ Ldω) ∼ 10^{2}). If the system size is smaller and/or the environment is more viscous, the overdamped approximation applies further strictly.

↵2 The absence of directional motion generally holds for any internal tensiondriven movement of a onedimensional line object on a wet substrate with spatially uniform η. Here, we prove this for a system of finite size with tensionfree ends. The equation of motion of such an object is generally written as ηu_{t} = T_{s}, where T(s,t) is the internal tension. The average velocity V of the system is given by , where L and R represent the values of s at the left and righthand boundaries, respectively. Because η is independent of s (but may depend on t in this argument), V ∝ (T(L) − T(R)) = 0 for free ends, fulfilling the condition that T(L) = T(R) = 0. We can arrive at the same consequence, V = 0, for the case of an infinite system size with a spatial period λ, because T(s) = T(s + λ) for any s.

↵3 The somewhat peculiar dimension of M, (m(J m^{−1})^{−1}), is a consequence of the choice of the averaged power density as the denominator of the definitional equation of M (equation (5.4)). The reason why we adopt the power density rather than the power is as follows: the power inputted into a region of the system inevitably increases with the length of the region, i.e. the interval λ of the above integration of (the longer the region, the more the energy dissipates). Thus, to make a fair comparison for the mileage quantity, it is reasonable to renormalize the power by the length of the region. We emphasize that after the normalization by 1/λ, the mileage M defined by equation (5.4) depends on the wavenumber k = 2π /λ, representing how many oscillations exist in a unit length of the system (figure 4c).

↵4 A small test is set here: can the propagation of the elongation–contraction wave cause locomotion in the case where the friction modulation contains only the random component, P(s,t) = R(s) or P(s,t) = R(X(s,t))?
 Received June 1, 2011.
 Accepted July 19, 2011.
 This journal is © 2011 The Royal Society