Generation of Direct-, Retrograde-, and Source-Wave Gaits in Multi-Legged Locomotion in a Decentralized Manner via Embodied Sensorimotor Interaction

Abstract

Multi-legged animals show several types of ipsilateral interlimb coordination. Millipedes use a direct-wave gait, in which the swing leg movements propagate from posterior to anterior. In contrast, centipedes use a retrograde-wave gait, in which the swing leg movements propagate from anterior to posterior. Interestingly, when millipedes walk in a specific way, both direct and retrograde waves of the swing leg movements appear with the waves' source, which we call the source-wave gait. However, the gait generation mechanism is still unclear because of the complex nature of the interaction between neural control and dynamic body systems. The present study used a simple model to understand the mechanism better, primarily how local sensory feedback affects multi-legged locomotion. The model comprises a multi-legged body and its locomotion control system using biologically inspired oscillators with local sensory feedback, phase resetting. Each oscillator controls each leg independently. Our simulation produced the above three types of animal gaits. These gaits are not predesigned but emerge through the interaction between the neural control and dynamic body systems through sensory feedback (embodied sensorimotor interaction) in a decentralized manner. The analytical description of these gaits' solution and stability clarifies the embodied sensorimotor interaction's functional roles in the interlimb coordination.

1. Introduction

Multi-legged animals, even those with a large number of legs, use several types of ipsilateral interlimb coordination according to the species and situations. Centipedes use a retrograde-wave gait in which the swing leg movements propagate from the anterior to posterior (Full, 1997; Kuroda et al., 2014). In contrast, millipedes generally use a direct-wave gait, in which the swing leg movements propagate from the posterior to anterior (Full, 1997; Kuroda et al., 2014). More interestingly, it is reported that when millipedes walk with the body axis bent like a U shape, both direct and retrograde waves of the swing leg movements appear at the source of the waves (Tamura et al., 2016), which we call the source-wave gait. However, it remains unclear what mechanisms generate these different types of interlimb coordination in multi-legged locomotion.

Regarding the locomotion of insects and mammals that use a direct-wave gait, where the swing leg movements propagate from the posterior to anterior, physiological studies have suggested that central pattern generators (CPGs) and sensory feedback play important roles in the interlimb coordination (Delcomyn, 2004; Büschges et al., 2008; Ijspeert, 2008). The CPGs generate rhythmic outputs, which are modulated by sensory feedback. The sensory feedback is critical especially during slow walking, as seen in stick insects (Delcomyn, 2004; Büschges et al., 2008).

However, the mechanism of the interlimb coordination has not been fully understood only from physiological studies because the locomotion is generated through complex interactions between motor control and body dynamics through sensory feedback (embodied sensorimotor interaction). Thus, many mathematical models and robots have been developed to clarify the functional roles of these interactions in the interlimb coordination (Steingrube et al., 2010; Owaki et al., 2012; Aoi et al., 2013; Schilling et al., 2013; Ambe et al., 2015).

Recently, mathematical models and robots for multi-legged locomotion have reproduced many aspects of the gaits of multi-legged animals by focusing on the embodied sensorimotor interaction. Tamura et al. (2016) showed that sensory feedback of load information generates a millipede-like gait. Yasui et al. (2017) and Kano et al. (2017) proposed distributed control schemes with sensory feedback of load information to reproduce the multi-legged locomotion observed when a part of the terrain is removed. They also investigated the transition between swimming and walking of centipedes based on the interplay of brain, CPG, and sensory feedback (Yasui et al., 2019). However, the mechanism for generating various interlimb coordination in multi-legged locomotion is still not fully understood.

The purpose of this study was to use a simple model to demonstrate that local sensory feedback generates the various interlimb coordination observed in multi-legged animals via embodied sensory-motor interactions. For that purpose, we constructed a mechanical model that imitates the flexible body of multi-legged animals and a control model that uses phase oscillators inspired by CPGs and phase resetting as local sensory feedback. The simulation results show that although the oscillators do not interact directly, three types of gaits (direct-, retrograde-, and source-wave gaits) emerge through the embodied sensorimotor interaction. Furthermore, we derive these analytical solutions and stabilities under some assumptions, which produce the three types of gaits regardless of the number of legs and clarify how the local sensory feedback generates these gaits through the embodied sensorimotor interaction.

2. Simulation

2.1. Mechanical Model

Skeletal structure and muscle arrangement have large effects on animal locomotion (Ting and Chiel, 2017), and appropriate models need to be created depending on the motion of interest. In particular, in fast locomotion where inertial effects are larger than viscous ones, body elasticity mainly determines the motion. For example, a spring-loaded inverted pendulum (SLIP) model has been widely used to investigate the characteristics of running in mammals (Full and Koditschek, 1999; Tanase et al., 2015; Adachi et al., 2020; Kamimura et al., 2021) and in cockroaches (Seipel et al., 2004; Spence et al., 2010). Detailed musculoskeletal models have been also used (Proctor and Holmes, 2018). Conversely, in slow locomotion, viscosity plays a dominant role. In stick insects, soon after swing muscle activity stops, swing leg movement ceases (Hooper et al., 2009). The effects of sensory feedback are dominant for slow-walking insects (Daun-Gruhn and Büschges, 2011).

In this study, we focus on relatively slow locomotion to investigate the functional roles of sensory feedback in multi-legged locomotion. We construct a simple mechanical model with flexible body and legs. Specifically, the mechanical model is composed of (N + 1) mass points, whose mass is m_i and whose height is x_i (i = 1, …, N + 1), as shown in Figure 1. The mass points move only vertically and are connected by springs and dampers (all the spring constants are κ and the damper constants are σ). The neutral length of the spring is 0. Each mass point has a massless leg (Leg i), which also moves vertically and whose length is l_i. The ground is modeled as a spring and damper (with spring constant K and damper constant D). The ground is much stiffer than the body spring (K ≫ κ) and the damper coefficient is set to provide overdamping. We used m_i = m for i = 2, …, N and m₁ = m_N+1 = m_B for the edges, where m/2 < m_B < m. The gravitational acceleration is g.

FIGURE 1

Figure 1. Mechanical model composed of (N+1) mass points. The mass points move only vertically and are connected via springs and dampers. Dotted lines represent the vertical positions of the mass points.

For N ≥ 3, the equations of motion are given by

mix¨i=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪−κ(x1−x2)−σ(x˙1−x˙2)−m1g −p1{K(x1−l1)+D(x˙1−l˙1)}−κ(2xi−xi+1−xi−1)−σ(2x˙i−x˙i+1−x˙i−1)−mig −pi{K(xi−l)+D(x˙i−l˙i)}−κ(xN+1−xN)−σ(x˙N+1−x˙N)−mN+1g −pN+1{K(xN+1−lN+1)+D(x˙N+1−l˙N+1)}i=1i=2,…,Ni=N+1 (1) mix¨i={−κ(x1−x2)−σ(x˙1−x˙2)−m1g −p1{K(x1−l1)+D(x˙1−l˙1)}i=1−κ(2xi−xi+1−xi−1)−σ(2x˙i−x˙i+1−x˙i−1)−mig −pi{K(xi−l)+D(x˙i−l˙i)}i=2,…,N−κ(xN+1−xN)−σ(x˙N+1−x˙N)−mN+1g −pN+1{K(xN+1−lN+1)+D(x˙N+1−l˙N+1)}i=N+1 (1)

where p_i(i = 1, …, N + 1) represents whether Leg i is in contact with the ground and is defined by,

pi={01li<xiotherwise. (2) pi={0li<xi1otherwise. (2)

2.2. Control Model

To understand the mechanism that generates rhythmic leg movement, central pattern generators (CPGs) have been well studied (Ijspeert, 2008; Daun-Gruhn and Büschges, 2011) and many modeling approaches are available for CPGs, such as relatively detailed biophysical models based on Hodgkin-Huxley-type neuron models (Rybak et al., 2006; Daun et al., 2009), connectionist models composed of simplified neuron models (e.g., leaky-integrator neurons) (Ijspeert, 2001; Pasemann et al., 2003) and abstract models using van der Pol and Matsuoka oscillators (van der Pol, 1926; Matsuoka, 1987). In contrast, some studies proposed reflex chains in place of CPGs to generate rhythmic leg movement, such as Walknet (Cruse et al., 1998), and others proposed heteroclinic oscillator models to represent intermediate behavior between the CPGs and reflex chains (Shaw et al., 2015).

These models are nonlinear, interact with sensory feedback, and show remarkable adaptation to the environment and prominent replication of animal movement, such as entrainment to limit cycle and extending the timing of the state transition. However, these models are complicated to understand the functional role of embodied sensorimotor interaction for interlimb coordination. We use a simple phase oscillator model developed in our previous studies (Aoi et al., 2017; Ambe et al., 2018) for better understanding of the functional role through the analytical description. Specifically, each leg has a phase oscillator whose phase is ϕ_i ∈ [0, 2π) (mod 2π) (i = 1, …, N+1). The leg length l_i is determined by ϕ_i as follows (Figure 2):

li=l(ϕi)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪LL−aϕi−2βπ(1−β)πL−a2π−ϕi(1−β)π0≤ϕi≤2βπ2βπ<ϕi<(1+β)π(1+β)π≤ϕi<2π (3) li=l(ϕi)={L0≤ϕi≤2βπL-aϕi-2βπ(1-β)π2βπ<ϕi<(1+β)πL-a2π-ϕi(1-β)π(1+β)π≤ϕi<2π (3)

where L, a(< L), and β ∈ (0.5, 1) are the maximum length, the amplitude of the swing leg movement, and the duty factor (ratio between the stance phase and step cycle durations), respectively. We defined the leg length by a piecewise-linear function and used mg/κ < a so that the leg is in the air in the swing phase (ϕ_i ∈ (2βπ, 2π)). While the length is L in the stance phase (ϕ_i ∈ [0, 2βπ]), it shortens in the takeoff phase (first half of the swing phase, ϕ_i ∈ (2βπ, (1+β)π)) and lengthens in the landing phase (second half of the swing phase, ϕ_i ∈ [(1 + β)π, 2π)). The oscillator phase follows the following dynamics:

dϕidt=ω+yi, (4) dϕidt=ω+yi, (4)

where ω is the basic frequency of locomotion and y_i represents a sensory feedback. We used ω≪κ/mi−−−−√ ω≪κ/mi to make the model walk slowly.

FIGURE 2

Figure 2. (A) Phase oscillator model. (B) Leg length li depending on phase [phi]i.

Sensory feedback has been physiologically well studied in insects (Delcomyn, 2004; Büschges et al., 2008). Insects have mechanoreceptors that can sense various information, such as contact with the ground, joint angles, and mechanical load at the nearby leg joint (Tuthill and Wilson, 2016). Especially in stick insects, sensory feedback plays an important role for inter-joint coordination. Strain signals from the trochanter play a major role in shaping thorax-coxa (TC)-joint motoneuronal activity during walking and contribute to the coordination of the TC-joint movement with the stepping pattern of the distal leg joints (Akay et al., 2004). Centipedes also sense load information and stop their periodic leg movements when the terrain is removed (Yasui et al., 2017). Locomotion rhythm resetting and phase shifting in motorneuron activities by sensory feedback and perturbation (phase resetting) has been observed in insects (Büschges, 1995), as well as in mammals (Schomburg et al., 1998; Rybak et al., 2006). Because phase resetting is amenable to mathematical analysis to determine the essential functional role of sensory feedback, this study focuses on phase resetting based on load information. Based on our previous studies (Aoi et al., 2013; Ambe et al., 2015, 2018), we incorporated the phase resetting mechanism as

τdyidt=−yi+ui (5) τdyidt=-yi+ui (5)

ui={0(2π−ϕi)δ(t−tii)0≤ϕi<(1+β)π(1+β)π≤ϕi<2π (6) ui={00≤ϕi<(1+β)π(2π-ϕi)δ(t-tii)(1+β)π≤ϕi<2π (6)

where tii tii is the time when Leg i touches the ground and δ() is Dirac's delta function. When Leg i touches the ground in the landing phase ((1 + β)π ≤ ϕ_i < 2π), the phase ϕ_i is reset to zero as shown in Figure 2A. In the present study, we used a first-order lag system with time constant τ to change the phase value continuously after the phase resetting for the simulation.

Because the leg movements of our model are determined by the oscillation phases, the relative phases between the oscillators ψ_i ∈ [0, 2π) (mod 2π)(i = 1, …, N) explain the gait, which is given by

ψi=ϕi+1−ϕi. (7) ψi=ϕi+1-ϕi. (7)

We investigated where the relative phases converged through the mechanical dynamics (1) and phase dynamics (4).

2.3. Results

We conducted forward dynamic simulations of our model in case of N = 3 to find stable gaits. All simulations used the parameters in Table 1. To investigate the characteristics of the gaits, we defined a Poincaré section Σ¹ for the relative phases just before Leg 1 touches the ground, where we used Ψ1≡[ψ11ψ12ψ13]T Ψ1≡[ψ11ψ21ψ31]T and ()ⁱ represents the value just before Leg i touches the ground.

TABLE 1

Table 1. Parameters for simulation.

2.3.1. Periodic Solutions

First, we used six sets of the initial relative phases for Ψ¹ to investigate if and how they converged. Figure 3 shows the time profile of the relative phases on Σ¹ for the six sets of the initial relative phases. Some of the relative phases converged to identical values, but others converged to different values. These converged relative phases correspond to periodic solutions. These results suggest that there are many periodic solutions.

FIGURE 3

Figure 3. Time profile of relative phases for Ψ1 simulated from six sets of initial values: (A) ψ11 ψ11 , (B) ψ12 ψ21 , and (C) ψ13 ψ31 .

To illustrate details of the periodic solutions, Figure 4A uses thick colored lines (red, blue, and green) to show the relative phases converged from many initial values (25 × 25 × 25) for Ψ¹. The solutions generally consist of three connected orthogonal line segments and lay close to two planes parallel to the ψ11 ψ11 – ψ12 ψ21 and ψ12 ψ21 – ψ13 ψ31 planes. Figures 4B,C show the solutions projected onto each plane. Near the connections between line segments, the segments are slightly bent. The three line segments correspond to the source-wave, retrograde-wave, and direct-wave gaits, as investigated below. The thin black line segments represent the analytical approximate solutions obtained in section 3.

FIGURE 4

Figure 4. Relative phases for Ψ1 of periodic solutions: (A) shown in ψ11 ψ11 - ψ12 ψ21 - ψ13 ψ31 space, (B) projected to a plane parallel to the ψ11 ψ11 - ψ12 ψ21 plane, and (C) projected to a plane parallel to the ψ12 ψ21 - ψ13 ψ31 plane. Thick colored lines and thin black lines represent simulated and analytic solutions, respectively. Points a, b, and c are used to show footprint diagrams in Figure 6.

Figures 5A–C show the slices of basin of attraction under the conditions ψ13=2(1−β)π,π ψ31=2(1-β)π,π , and 2βπ, respectively. All the initial values (99 × 99) converged to one of the three gaits, and the basin of attraction was separated into three parts, each corresponding to one gait type.

FIGURE 5

Figure 5. Slices of basin of attraction for three gaits in simulation under three conditions: (A) ψ13=2(1−[beta])[pi] ψ31=2(1-[beta])[pi] , (B) ψ13=[pi] ψ31=[pi] , and (C) ψ13=2[beta][pi] ψ31=2[beta][pi] .

2.3.2. Characteristics of Three Types of Gaits

The obtained periodic solutions were categorized by three gaits: direct-, retrograde-, and source-wave gaits. The direct-wave gait has the following relative phases:

ψˆ11≈2(1−β)π,ψˆ12≈2(1−β)π,2(1−β)π⪅ψˆ13⪅2βπ ψ^11≈2(1-β)π,ψ^21≈2(1-β)π,2(1-β)π⪅ψ^31⪅2βπ

where ()ˆ ()^ indicates the periodic solution. This gait appears as one line segment, as shown in Figure 4. Figure 6A shows the footprint diagram for Ψˆ1=[1.13 0.98 1.50]T Ψ^1=[1.13 0.98 1.50]T (point a in Figure 4). In this gait, the swing leg movement of Legs 1, 2, and 3 propagates from posterior to anterior with the same interval, as shown by red arrows. However, the relative phase of the swing leg movement between Legs 3 and 4 can vary within the range indicated by the violet arrow.

FIGURE 6

Figure 6. Footprint diagrams of periodic solutions: (A) Direct-wave gait (point a in Figure 4), (B) retrograde-wave gait (point b in Figure 4), and (C) source-wave gait (point c in Figure 4). Black bars can change to gray bars as long as the violet box is in the range indicated by the violet arrow.

The retrograde-wave gait has the following relative phases:

2(1−β)π⪅ψˆ11⪅2βπ,ψˆ12≈2βπ,ψˆ13≈2βπ. 2(1-β)π⪅ψ^11⪅2βπ,ψ^21≈2βπ,ψ^31≈2βπ.

This also appears as one line segment (to be precise, it is slightly curved) in Figure 4. Figure 6B shows the footprint diagram for Ψˆ1=[3.97 5.29 5.29]T Ψ^1=[3.97 5.29 5.29]T (point b in Figure 4). In this gait, the swing leg movement of Legs 2, 3, and 4 propagates from anterior to posterior with the same interval, as shown by green arrows. However, the relative phase of the swing leg movement between Legs 1 and 2 can vary in the range indicated by the violet arrow.

The source-wave gait has the following relative phases:

ψˆ11≈2(1−β)π,2(1−β)π⪅ψˆ12⪅2βπ,ψˆ13≈2βπ. ψ^11≈2(1-β)π,2(1-β)π⪅ψ^21⪅2βπ,ψ^31≈2βπ.

This also appears as one line segment (to be precise, it comprises two line segments) in Figure 4. Figure 6C shows the footprint diagram for Ψ¹ = [1.13 1.10 5.28]^T (point c in Figure 4). In this gait, the swing leg movement of Legs 1 and 2 propagates from posterior to anterior, as shown by red arrows, while that of Legs 3 and 4 propagates from anterior to posterior as shown by green arrows. However, the relative phase of the swing leg movement between Legs 2 and 3 can vary in the range indicated by the violet arrow.

These three types of gaits were not predetermined, but emerged through the mechanical and sensory interaction. These gaits are represented by line segments for Ψ¹, which are serially connected as shown in Figure 4.

3. Analysis

Three types of gaits were obtained via the mechanical and sensory interaction in simulations with a multi-mass spring model. Although the models described in section 2 used four mass points, these gaits are expected to appear in models with more mass points (attached Supplementary Video 1). Figure 7 shows the three types of gaits for N+1 mass points. For the direct-wave gait (Figure 7A), the swing leg movements propagate anteriorly, where there is only one swing leg within one wavelength. While the phase difference between the neighboring legs for Legs 1 to N is around 2(1−β)π, that between Legs N and N+1 is not unique and is variable from 2(1−β)π to 2βπ. For the retrograde-wave gait (Figure 7B), the swing leg movements propagate posteriorly, where there is only one swing leg within one wavelength. While the phase difference between the neighboring legs for Legs 2 to N + 1 is around 2βπ, that between Legs 1 and 2 is not unique and is variable from 2(1 − β)π to 2βπ. For the source-wave gait (Figure 7C), while the anterior swing leg movements (Legs 1 to k − 1, k ∈ [2, N]) propagate anteriorly, those of the posterior part (Legs k to N + 1) propagate posteriorly. Although the phase differences of the neighboring legs for the anterior and posterior parts are around 2(1 − β)π and 2βπ, respectively, that between Legs k − 1 and k is not unique and is variable from 2(1 − β)π to 2βπ. To clarify the mechanism for generating these three gaits, we simplified our model by some physical assumptions to derive analytical solutions.

FIGURE 7

Figure 7. Three types of gaits for N+1 mass points: (A) Direct-wave, (B) retrograde-wave, and (C) source-wave gaits. In the direct-wave (retrograde-wave) gait, while the swing leg movements of Legs 1 to N (Legs 2 to N+1) propagate anteriorly (posteriorly) with a constant relative phase, the relative phase between Legs N and N+1 (Legs 1 and 2) is not unique. In the source-wave gait, while the swing leg movements of the anterior part (Legs 1 to k−1, k [element of] [2, N]) propagate anteriorly, those of the posterior part (Legs k to N + 1) propagate posteriorly. The relative phase between Legs k−1 and k is not unique. In these three gaits, there is only one swing leg within one wavelength.

A concrete way to derive the analytical solutions is as follows. First, we simplify our mechanical and controller models. In particular, we assume that the natural frequencies of the mass points are larger than the gait frequency so that the vertical position of the mass points x_i (i = 1, …, N + 1) is determined uniquely by the oscillator phase ϕ_i in section 3.1. We also assume that the time constant of our control model is smaller than the gait cycle duration found in section 3.2, which allows the Poincaré map to be obtained by determining the phase resetting value ϕii ϕii . These assumptions give a touchdown relation between ϕii,ϕii−1 ϕii,ϕi-1i , and ϕii+1 ϕi+1i . In section 3.3, we determine the order of the touchdown events in the case of the source-wave gait based on the simulation results and rewrite the touchdown relation by ϕi−1i−1,ϕii,ϕi+1i+1 ϕi-1i-1,ϕii,ϕi+1i+1 , and the relative phases between the oscillators, which gives periodic solutions of the source-wave gait. In section 3.4, we investigate the stability of the periodic solutions. In section 3.5, we derive the periodic solutions and stability of the direct-wave gait from the front–rear symmetry. In section 3.6, we compare the analytical results with the simulation results. In section 3.7, we describe the geometrical meaning of the obtained solutions.

3.1. Simplification of Mechanical Model

We used ω≪κ/mi−−−−√ ω≪κ/mi in (4) to make the model walk slowly in the simulation. Thus, we ignored the dynamics of the mass points and focused only on the equilibrium of forces by replacing the ground reaction forces with unknown variables. Then, (1) becomes

⎧⎩⎨⎪⎪κ(x2−x1)+R1=m1gκ(xi−1−xi)+κ(xi+1−xi)+Ri=migκ(xN−xN+1)+RN+1=mN+1gi=1i∈[2,N]i=N+1, (8) {κ(x2-x1)+R1=m1gi=1κ(xi-1-xi)+κ(xi+1-xi)+Ri=migi∈[2,N]κ(xN-xN+1)+RN+1=mN+1gi=N+1, (8)

where R_i is the ground reaction force. The number of unknown variables is 2(N + 1); x₁, …, x_N+1 and R₁, …, R_N+1. These variables depend on the foot contact conditions. In particular, when Leg i is in the air, x_i ≥ l(ϕ_i) and R_i = 0. In contrast, when Leg i is in contact with the ground, x_i = l(ϕ_i) and R_i > 0. Therefore, when ϕ_i and foot contact conditions for all legs are given, the number of the unknown variables is reduced to N + 1 and the unknown variables are solved by (8).

3.2. Simplification of Control Model

Although we used a first-order lag system in (9) for the simulation, the time constant τ is small compared with the gait period. Therefore, we set τ = 0, which reduces (4) to

dϕidt=ω+ui. (9) dϕidt=ω+ui. (9)

This gives

dψidt=ui+1−uii=1,…,N. (10) dψidt=ui+1-uii=1,…,N. (10)

A Poincaré section Σ^j (j = 1, …, N+1) was defined for the relative phases just before the touchdown of Leg j and we used Ψj≡[ψj1ψj2…ψjN]T Ψj≡[ψ1jψ2j…ψNj]T . The Poincaré map for Σ^j was given by integrating (10) for one period. Under the assumption that each leg experiences phase resetting once each period, the map is given by

ψji↦ψji+ϕii−ϕi+1i+1i=1,…,N. (11) ψij↦ψij+ϕii-ϕi+1i+1i=1,…,N. (11)

To obtain the Poincaré map, ϕii(i=1,…,N+1) ϕii(i=1,…,N+1) needs to be determined.

Based on the periodic solutions, we assume that there is only one swing leg within one wavelength, just before Leg i touches the ground (i = 1, …, N+1), the neighboring legs (Legs i−1 and i+1) are in contact with the ground. Therefore, xi−1=l(ϕii−1) xi-1=l(ϕi-1i) , xi=l(ϕii) xi=l(ϕii) , xi+1=l(ϕii+1) xi+1=l(ϕi+1i) , R_i−1 > 0, R_i = 0, and R_i+1 > 0 are satisfied just before the touchdown of Leg i. By substituting these variables into (8), we obtain

l(ϕii)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪fB(ϕi2)f(ϕii−1,ϕii+1)fB(ϕiN)i=1i∈[2,N]i=N+1, (12) l(ϕii)={fB(ϕ2i)i=1f(ϕi-1i,ϕi+1i)i∈[2,N]fB(ϕNi)i=N+1, (12)

where

f(ϕh,ϕf)=12{l(ϕh)+l(ϕf)}−mg2κ (13) f(ϕh,ϕf)=12{l(ϕh)+l(ϕf)}-mg2κ (13)

fB(ϕ)=l(ϕ)−mBgκ. fB(ϕ)=l(ϕ)-mBgκ.

The substitution of (12) into (3) gives one relation between ϕii−1,ϕii ϕi-1i,ϕii , and ϕii+1 ϕi+1i for each i (i = 1, …, N + 1). Because we can write ϕii−1 ϕi-1i and ϕii+1 ϕi+1i as ϕi−1i−1,ϕii,ϕi+1i+1 ϕi-1i-1,ϕii,ϕi+1i+1 , and Ψ^j, as discussed in the following section, ϕii ϕii is obtained by Ψ^j. As a result, we obtain the Poincaré map from (11).

3.3. Analytical Solution of Source-Wave Gait

Here, we derive the periodic solution of the source-wave gait characterized by Leg k ∈ [2, N] that determines the boundary between the direct- and retrograde-wave regions (Figure 7C). The solution of the source-wave gait is useful for deriving those of the direct- and retrograde-wave gaits. In particular, the solution of the retrograde-wave gait is obtained by k = 2 for the solution of the source-wave gait (Figure 7B). However, note that k = N+1 does not fully explain the solution of the direct-wave gait, which is instead derived by the front–rear symmetry and the solution of the retrograde-wave gait, as explained in section 3.5.

To derive the periodic solution of the source-wave gait, we first define the Poincaré section as Σ^k and transform the right-hand side of (12) to be represented by the relative phases Ψ^k and the touchdown phase ϕii ϕii (i = 1, …, N + 1). For that purpose, we assume some conditions for the phases of the neighboring legs at each touchdown event, determine the order of the touchdown events, and represent the relative phases of the neighboring legs of the touchdown leg using Ψ^k and ϕii ϕii by accounting for the effects of phase resetting. Second, we derive the conditions of ϕii ϕii for the periodicity. The reduced equation of (12) and the periodicity conditions yield the periodic solution.

Based on the simulation results obtained in section 2, we assume the following conditions for the phases of the neighboring legs at each touchdown event (Figure 8). For the touchdown of Leg k (boundary between direct- and retrograde-wave regions), ϕkk−1 ϕk-1k is in the stance or takeoff phase ( 0≤ϕkk−1<(1+β)π 0≤ϕk-1k<(1+β)π ), and ϕkk+1 ϕk+1k is in the takeoff phase ( 2βπ≤ϕkk+1<(1+β)π 2βπ≤ϕk+1k<(1+β)π ). For the touchdown of Leg j ∈ [2, k − 1] (direct-wave region), ϕjj−1 ϕj-1j is in the takeoff phase ( 2βπ≤ϕjj−1<(1+β)π 2βπ≤ϕj-1j<(1+β)π ) and ϕjj+1 ϕj+1j is in the stance phase ( 0≤ϕjj+1<2βπ 0≤ϕj+1j<2βπ ) because the swing movements of Legs 1 to k − 1 propagate anteriorly. For the touchdown of Leg i ∈ [k + 1, N] (retrograde-wave region), ϕii−1 ϕi-1i is in the stance phase ( 0≤ϕii−1<2βπ 0≤ϕi-1i<2βπ ) and ϕii+1 ϕi+1i is in the takeoff phase ( 2βπ≤ϕii+1<(1+β)π 2βπ≤ϕi+1i<(1+β)π ) because the swing movements of Legs k to N + 1 propagate posteriorly. For the touchdown of Legs 1 and N + 1 (edges of direct- and retrograde-wave regions), ϕ12 ϕ21 and ϕN+1N ϕNN+1 are in the stance phases ( 0≤ϕ12,ϕN+1N<2βπ 0≤ϕ21,ϕNN+1<2βπ ).

FIGURE 8

Figure 8. Assumptions for the phases of the neighboring legs at each touchdown event for the source-wave gait.

These assumptions determine the order of touchdown events. In the retrograde-wave region, suppose that n(> k) is the minimum value that satisfies ϕkn∈(ϕn−1n,ϕnn] ϕnk∈(ϕnn-1,ϕnn] , which means that Legs (n−1), k, and n touch the ground sequentially (Figure 9A). Legs k to n correspond to almost one wavelength of the gait. Because the swing movement propagates posteriorly from Legs k to N−1, Leg n is in the swing phase when Leg k touches down. Thus, the order of touchdown events for Legs k to n is determined as being in the order k, n, k−1, k−2, k−3, … , n−1, as shown in Figure 9A. The order of touchdown events in the direct-wave region is determined similarly.

FIGURE 9

Figure 9. Evolution of oscillator phases within one wavelength: (A) Order of touchdown events of Legs k to n in the retrograde-wave region. Legs k, n, k+1, k+2, k+3, …, n−1 touch down sequentially. (B) Evolution of oscillator phases and relative phases at touchdown events of Legs k to n. Phase resetting at each touchdown changes only the oscillator phase of the touchdown leg and the relative phases of the neighboring legs. ψii−1 ψi-1i and ψii ψii just before the touchdown of Leg i (i = k, n, k+1, k+2, …, n−2, n−1) are represented by Ψk, [phi]i−1i−1 [phi]i-1i-1 , and [phi]i+1i+1 [phi]i+1i+1 , as highlighted in orange.

When the order of touchdown events is determined, the relative phases ψii−1 ψi-1i and ψii ψii just before touchdown of Leg i can be represented using Ψ^k and ϕi−1i−1 ϕi-1i-1 , and ϕi+1i+1 ϕi+1i+1 by accounting for the effects of phase resetting. Specifically, because phase resetting at each touchdown changes only the oscillator phase of the touchdown leg and the relative phases of the neighboring legs, ψii−1 ψi-1i and ψii ψii just before touchdown of Leg i can be represented using Ψ^k, ϕi−1i−1 ϕi-1i-1 , and ϕi+1i+1 ϕi+1i+1 , depending on the experiences of the leg touchdowns. In the case of the retrograde-wave region, the evolution of the relative phases for Legs k to n at each touchdown event can be written as Figure 9B. For i = k and n, ψii−1 ψi-1i and ψii ψii are given by ψii−1=ψki−1 ψi-1i=ψi-1k and ψii=ψki ψii=ψik , respectively (highlighted in orange at the touchdowns of Legs k and n in Figure 9B), because both neighboring legs do not experience phase resetting in the period between the touchdown events of Legs k and n. As a result, ϕii−1=ϕii−ψki−1 ϕi-1i=ϕii-ψi-1k and ϕii+1=ϕii+ψki−2π ϕi+1i=ϕii+ψik-2π ( ϕii+1 ϕi+1i is subtracted by 2π to satisfy 0≤ϕii+1<2π 0≤ϕi+1i<2π ) are satisfied because ψii−1=ϕii−ϕii−1 ψi-1i=ϕii-ϕi-1i and ψii=ϕii+1−ϕii ψii=ϕi+1i-ϕii . For i ∈ [k+1, n−2], ψii−1 ψi-1i and ψii ψii are given by ψii−1=ψki−1−(2π−ϕi−1i−1) ψi-1i=ψi-1k-(2π-ϕi-1i-1) and ψii=ψki ψii=ψik , respectively (highlighted in orange at the touchdown of Leg k+1 for i = k+1 in Figure 9B), because only the posterior Leg i−1 experiences phase resetting in the period between the touchdown events of Legs k and i. As a result, ϕii−1=ϕii−ψki−1+2π−ϕi−1i−1 ϕi-1i=ϕii-ψi-1k+2π-ϕi-1i-1 and ϕii+1=ϕii+ψki−2π ϕi+1i=ϕii+ψik-2π are satisfied. For i = n−1, ψii−1 ψi-1i and ψii ψii are given by ψii−1=ψki−1−(2π−ϕi−1i−1) ψi-1i=ψi-1k-(2π-ϕi-1i-1) and ψii=ψki+(2π−ϕi+1i+1) ψii=ψik+(2π-ϕi+1i+1) , respectively (highlighted in orange at the touchdown of Leg n−1 in Figure 9B), because both neighboring legs experience phase resetting in the interval between the touchdown events of Legs k and n−1. As a result, ϕii−1=ϕii−ψki−1+2π−ϕi−1i−1 ϕi-1i=ϕii-ψi-1k+2π-ϕi-1i-1 and ϕii+1=ϕii+ψki−ϕi+1i+1 ϕi+1i=ϕii+ψik-ϕi+1i+1 are satisfied. This means that ϕii−1 ϕi-1i and ϕii+1 ϕi+1i can be represented by Ψ^k, ϕi−1i−1 ϕi-1i-1 , ϕii ϕii , and ϕi+1i+1 ϕi+1i+1 . These analyses are also applicable to the direct-wave region.

By using these results, the right-hand side of (12) can be rewritten using ψki(i=1,…,N) ψik(i=1,…,N) and ϕii(i=1,…,N+1) ϕii(i=1,…,N+1) as

l(ϕii)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪fB(ϕ11+ψk1−2π)fB(ϕ11+ψk1−ϕ22)f(ϕii−ψki−1,ϕii+ψki−2π)f(ϕii−ψki−1+2π−ϕi−1i−1,ϕii+ψki−ϕi+1i+1)f(ϕii−ψki−1,ϕii+ψki−ϕi+1i+1)f(ϕii−ψki−1+2π−ϕi−1i−1,ϕii+ψki−2π)fB(ϕN+1N+1−ψkN)fB(ϕN+1N+1−ψkN+2π−ϕNN)i∈Sk1i∈Sk2i∈Sk3i∈Sk4i∈Sk5i∈Sk6i∈Sk7i∈Sk8 (14) l(ϕii)={fB(ϕ11+ψ1k-2π)i∈S1kfB(ϕ11+ψ1k-ϕ22)i∈S2kf(ϕii-ψi-1k,ϕii+ψik-2π)i∈S3kf(ϕii-ψi-1k+2π-ϕi-1i-1,ϕii+ψik-ϕi+1i+1)i∈S4kf(ϕii-ψi-1k,ϕii+ψik-ϕi+1i+1)i∈S5kf(ϕii-ψi-1k+2π-ϕi-1i-1,ϕii+ψik-2π)i∈S6kfB(ϕN+1N+1-ψNk)i∈S7kfB(ϕN+1N+1-ψNk+2π-ϕNN)i∈S8k (14)

where Sk1−8 S1-8k are classified based on whether the neighboring legs (Legs i − 1 and i + 1) of Leg i experience phase resetting in the interval between the touchdown events of Legs k and i, as shown in Table 2. Specifically, Sk1,2,7,8 S1,2,7,8k are for i = 1 and N + 1, which are given by

Sk1={i∣i=1}∩{i∣ϕki∈(ϕi+1i,ϕii]}Sk2={i∣i=1}∩{i∣i∉Sk1}Sk7={i∣i=N+1}∩{i∣ϕki∈(ϕi−1i,ϕii]}Sk8={i∣i=N+1}∩{i∣i∉Sk7}, S1k={i∣i=1}∩{i∣ϕik∈(ϕii+1,ϕii]}S2k={i∣i=1}∩{i∣i∉S1k}S7k={i∣i=N+1}∩{i∣ϕik∈(ϕii-1,ϕii]}S8k={i∣i=N+1}∩{i∣i∉S7k},

Sk3,4,5,6 S3,4,5,6k are for i = 2 to N, which are given by

Sk3={i∣i=k}∪{i∈[2,k−1]∣ϕki∈(ϕi+1i,ϕii]} ∪{i∈[k+1,N]∣ϕki∈(ϕi−1i,ϕii]}Sk4={i∈[2,k−1]∣ϕki−1∈(ϕii−1,ϕi−1i−1]} ∪{i∈[k+1,N]∣ϕki+1∈(ϕii+1,ϕi+1i+1]}Sk5={i∈[2,k−1]∣i∉Sk3 and i∉Sk4}Sk6={i∈[k+1,N]∣i∉Sk3 and i∉Sk4}. S3k={i∣i=k}∪{i∈[2,k-1]∣ϕik∈(ϕii+1,ϕii]} ∪{i∈[k+1,N]∣ϕik∈(ϕii-1,ϕii]}S4k={i∈[2,k-1]∣ϕi-1k∈(ϕi-1i,ϕi-1i-1]} ∪{i∈[k+1,N]∣ϕi+1k∈(ϕi+1i,ϕi+1i+1]}S5k={i∈[2,k-1]∣i∉S3k and i∉S4k}S6k={i∈[k+1,N]∣i∉S3k and i∉S4k}.

From (11), periodic solutions must satisfy

ϕˆ11=ϕˆ22=⋯=ϕˆN+1N+1≡ϕˆtd. (15) ϕ^11=ϕ^22=⋯=ϕ^N+1N+1≡ϕ^td. (15)

This means that the touchdown phase of all oscillators equals ϕˆtd ϕ^td . From (14) and (15), we obtain the periodic solution by ϕˆtd ϕ^td and ψˆki(i=1,…,N) ψ^ik(i=1,…,N) as follows:

ϕˆtd=2π−mBb (16) ϕ^td=2π-mBb (16)

ψˆki=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪ψˆkk−12(1−β)π−(2mB−m)b2(1−β)π−(3mB−m)bψˆkk2βπ+(2mB−m)b2βπ+(3mB−m)bi=k−1i∈Tk1i∈[1,k−2] and i∉Tk1i=ki+1∈Tk2otherwise, (17) ψ^ik={ψ^k-1ki=k-12(1-β)π-(2mB-m)bi∈T1k2(1-β)π-(3mB-m)bi∈[1,k-2] and i∉T1kψ^kki=k2βπ+(2mB-m)bi+1∈T2k2βπ+(3mB-m)botherwise, (17)

where b = ((1−β)gπ)/(κa). The first row, the following two rows, and the remaining rows of the right-hand side of (17) represent the relative phases on the boundary, in the direct-wave region, and in the retrograde-wave region, respectively. The sets Tk1,2 T1,2k are given by

Tk1={i∈[1,k−2]|2βπ+(2mB−m)b<ϕˆki≤ϕˆtd}Tk2={i∈[k+2,N+1]|2βπ+(2mB−m)b<ϕˆki≤ϕˆtd}. (18) T1k={i∈[1,k-2]|2βπ+(2mB-m)b<ϕ^ik≤ϕ^td}T2k={i∈[k+2,N+1]|2βπ+(2mB-m)b<ϕ^ik≤ϕ^td}. (18)

Tk1 T1k is the set of legs (Leg i) in the direct-wave region whose neighboring legs (Legs i−1 and i+1) do not experience phase resetting between the touchdowns of Legs k and i (i.e., Sk3 S3k in the direct-wave region). Tk2 T2k is the set of legs (Leg i) in the retrograde-wave region whose neighboring legs (Legs i−1 and i+1) do not experience phase resetting between the touchdowns of Legs k and i (i.e., Sk3 S3k in the retrograde-wave region). ψˆkk−1 ψ^k-1k is not determined uniquely but satisfies

2(1−β)π−(3mB−m)b<ψˆkk−1≤2βπ. (19) 2(1-β)π-(3mB-m)b<ψ^k-1k≤2βπ. (19)

This non-uniqueness is because the length of Leg i is constant in the stance phase (i.e., l(ϕ_i) = L, which does not determine ϕ_i uniquely). ψˆkk ψ^kk is represented using ψˆkk−1 ψ^k-1k as

ψˆkk=⎧⎩⎨⎪⎪ψˆkk−1+(4β−2)π+(4mB−m)b2βπ+(3mB−m)bψA<ψˆkk−1<ψBψB≤ψˆkk−1≤2βπ (20) ψ^kk={ψ^k-1k+(4β-2)π+(4mB-m)bψA<ψ^k-1k<ψB2βπ+(3mB-m)bψB≤ψ^k-1k≤2βπ (20)

where ψ_A = 2(1 − β)π − (3m_B − m)b and ψ_B = 2(1 − β)π − m_Bb. Therefore, the relationship of ψˆkk−1 ψ^k-1k and ψˆkk ψ^kk is explained by two connected segments, as shown in Figure 10. We call these segments solution sets A and B. Note that solution set A is smaller than solution set B, as shown in Figure 4, because κ ≫ 1.

TABLE 2

Table 2. Classification of sets Sk1−8 S1-8k .

FIGURE 10

Figure 10. Relation between ψˆkk−1 ψ^k-1k and ψˆkk ψ^kk .

3.4. Stability of Source-Wave Gait

To discuss stability, we added perturbation Δψki(i=1,…,N) Δψik(i=1,…,N) to the obtained periodic solutions on Σ^k and evaluated the linear stability by analytically calculating the linear map of perturbation:

ΔΨk↦PkΔΨk, (21) ΔΨk↦PkΔΨk, (21)

where ΔΨk=[Δψk1Δψk2…ΔψkN]T ΔΨk=[Δψ1kΔψ2k…ΔψNk]T , and Pk∈RN×N Pk∈ℝN×N is the Jacobian matrix of the Poincaré map. We assumed that the perturbation is too small to change the order of touchdown events and the sets S_i.

When we write the perturbed touchdown phase of Leg i as ϕˆtd+Δϕii ϕ^td+Δϕii and substitute ψki=ψˆki+Δψki ψik=ψ^ik+Δψik and ϕii=ϕˆtd+Δϕii ϕii=ϕ^td+Δϕii into (14), Δϕii Δϕii (i = 1, …, N + 1) can be represented with Δψki Δψik . The matrix P_k is derived by substituting ϕii=ϕˆtd+Δϕii ϕii=ϕ^td+Δϕii into (11). The eigenvalues λ_Sk of P_k for the source-wave gait are derived depending on solution sets A and B by

λSk=⎧⎩⎨⎪⎪{1,1/2,2/3,…,2/3}{1,2/3,2/3,…,2/3}ψA<ψˆkk−1<ψBψB≤ψˆkk−1≤2βπ (22) λSk={{1,1/2,2/3,…,2/3}ψA<ψ^k-1k<ψB{1,2/3,2/3,…,2/3}ψB≤ψ^k-1k≤2βπ (22)

There is only one eigenvalue of 1 for both solution sets due to the non-uniqueness. However, the other eigenvalues are <1, which means that any initial points near the solution set will converge to the solution set.

3.5. Analytical Solution and Stability of Direct-Wave Gait

The solution of the direct-wave gait cannot be derived by substituting k = N + 1 for those of the source-wave gaits. This is because part of the assumption regarding the phases of neighboring legs (specifically, that ϕkk−1 ϕk-1k is in the takeoff phase, as in the middle figure of Figure 8) is not correct when k = N + 1 for the direct-wave gait. Thus, we derive the solution of the direct-wave gait from the front–rear symmetry and the solution of the retrograde-wave gait.

We denote the flow of the oscillator phase i (i = 1, …, N + 1) with the initial value (ϕ˜k1,ϕ˜k2,…,ϕ˜kN+1) (ϕ~1k,ϕ~2k,…,ϕ~N+1k) just before the touchdown of Leg k by Φi(t;ϕ˜k1,ϕ˜k2,…,ϕ˜kN+1) Φi(t;ϕ~1k,ϕ~2k,…,ϕ~N+1k) . From the front–rear symmetry of our model, the following equation is satisfied for the initial phases such that ϕ˜ki=ϕ˜N+2−ki ϕ~ik=ϕ~iN+2-k for i = 1, …, N + 1:

Φi(t;ϕ˜k1,ϕ˜k2,…,ϕ˜kN+1) =ΦN+2−i(t;ϕ˜N+2−kN+1,ϕ˜N+2−kN,…,ϕ˜N+2−k1) (23) Φi(t;ϕ~1k,ϕ~2k,…,ϕ~N+1k) =ΦN+2-i(t;ϕ~N+1N+2-k,ϕ~NN+2-k,…,ϕ~1N+2-k) (23)

The periodic solution of the direct-wave gait ψˆNi(i=1,…,N) ψ^iN(i=1,…,N) on Σ^N, where ϕˆtd ϕ^td is the same as (16), is derived from solution set B of the retrograde-wave gait (k = 2 for the source-wave gait) using the symmetry condition (symmetrical solution for solution set A of the retrograde-wave gait equals solution set A of the source-wave gait with k = N). Specifically, the substitution of k = N into (23) yields Φi(t;ϕ˜NN+1,…,ϕ˜N1)=ΦN+2−i(t;ϕ˜21,…,ϕ˜2N+1) Φi(t;ϕ~N+1N,…,ϕ~1N)=ΦN+2-i(t;ϕ~12,…,ϕ~N+12) . When the periodic solution ϕˆNj ϕ^jN (j = 1, …, N + 1) is used for the initial value ϕ˜Nj ϕ~jN , we obtain ϕˆNi=ϕˆ2N+2−i ϕ^iN=ϕ^N+2-i2 (i = 1, …, N). This yields ψˆNi=2π−ψˆ2N+1−i ψ^iN=2π-ψ^N+1-i2 because ψˆNi=ϕˆNi+1−ϕˆNi ψ^iN=ϕ^i+1N-ϕ^iN . By substituting solution set B of the retrograde-wave gait ψˆ2N+1−i ψ^N+1-i2 (k = 2 in (17)), we obtain the solution of the direct-wave gait as

ψˆNi=⎧⎩⎨⎪⎪⎪⎪ψˆNN2(1−β)π−(2mB−m)b2(1−β)π−(3mB−m)bi=Ni∈TN1otherwise. (24) ψ^iN={ψ^NNi=N2(1-β)π-(2mB-m)bi∈T1N2(1-β)π-(3mB-m)botherwise. (24)

where ψˆNN ψ^NN is an arbitrary constant fulfilling

2(1−β)π≤ψˆNN≤2βπ+mBb. (25) 2(1-β)π≤ψ^NN≤2βπ+mBb. (25)

The linear stability of the solution can also be calculated using the symmetry condition (23). The symmetry condition gives the relation ΔψNi=−Δψ2N+1−i ΔψiN=-ΔψN+1-i2 , which is represented as ΔΨ^N = −QΔΨ², where the matrix Q is an anti-diagonal matrix with all the anti-diagonal elements of 1. By using Q, the Jacobian matrix of the Poincaré map of the direct-wave gait P_D is calculated as PD=QP2Q−1 PD=QP2Q-1 . Thus, the eigenvalues λ_D of P_D are obtained by those of P₂ as

λD={1,2/3,2/3,…,2/3}. (26) λD={1,2/3,2/3,…,2/3}. (26)

There is only one eigenvalue of 1 due to the non-uniqueness. However, the other eigenvalues are <1, which means that any initial points near the solution set will converge to the solution set.

All the solution sets we derived were connected serially as a chain. Specifically, the boundary on solution set A of the source-wave gait with k = i (left side of Figure 10) and the boundary on solution set B of the source-wave gait with k = i + 1 (right side of Figure 10) are connected, as shown in Figure 11A. In addition, the boundary on solution set A of the source-wave gait with k = i + 1 and the boundary on solution set B of the source-wave gait with k = i + 2 are connected, as shown in Figure 11A. Furthermore, the boundary on solution set A of the source-wave gait with k = N and the boundary of the direct-wave gait are connected, as shown in Figure 11B. The boundary on solution set B of the source-wave gait with k = 3 and the boundary on solution set A of the source-wave gait with k = 2 (retrograde-wave gait) are connected, as shown in Figure 11C. Thus, the obtained solution sets consist of many connected segments constrained in different planes, and the boundaries of the whole solution set correspond to the direct- and retrograde-wave gaits. Note that discontinuous jumps may exist on the analytical solution sets, as highlighted in Figure 11A. The jumps occur when the elements of the sets T₁ and T₂ change, that is, when the order of touchdown events of the neighboring legs changes. For example, when Leg s becomes an element of Tk1 T1k , that is, when the order of touchdown events of the neighboring legs changes from Legs s+1, s, and s−1 to Legs s, s−1, and s+1, the relative phase ψˆks ψ^sk jumps from 2(1−β)π−(3m_B−m)b to 2(1−β)π−(2m_B−m)b, as shown in (17).

FIGURE 11

Figure 11. Connection of the solution sets: (A) Source-wave gaits (k = i, i + 1, i + 2) projected onto the ψˆii−1 ψ^i-1i - ψˆii ψ^ii - ψˆii+1 ψ^i+1i space; (B) direct- and source-wave gaits (k = N−1, N) projected onto the ψˆN−1N−2 ψ^N-2N-1 - ψˆN−1N−1 ψ^N-1N-1 - ψˆN−1N ψ^NN-1 space; and (C) retrograde- (k = 2) and source-wave gaits (k = 3, 4) projected onto the ψˆ21 ψ^12 - ψˆ22 ψ^22 - ψˆ23 ψ^32 space.

3.6. Comparison With Simulation Results

To validate the analytical results, we compared the obtained solutions for N = 3 with the simulation results on Σ¹ in Figure 4, where the analytical solutions were derived using the same parameters of the simulation. In both the simulation and analytical results, source- and retrograde-wave gaits show two segments as solution sets A and B (although some discontinuous jumps appear). The segments of three gaits are serially connected, and their edges correspond to the direct- and retrograde-wave gaits. These characteristics were identical between the analytical and simulation results. Furthermore, the maximum distance between the solution sets of the simulation and analytical results was only 0.11 (1.8% of 2π). The errors mainly came from discontinuous jumps in the analytical solution, which are not shown in the simulation, because the first-order lag system in (5) smoothed the jumps.

3.7. Geometrical Features of Solutions

From the analytical solution (17), we found that the swing movements of the source-wave gait are isotropic waves and propagate forward from Leg k − 1 (direct-wave region) and backward from Leg k (retrograde-wave region) with constant speeds. This gait is similar to that observed in millipedes (Tamura et al., 2016). When the interval between the legs is η, the wavelength ξ and velocity v of the swing movement are derived by

ξ=11−β−mBg2κa1−(3mB−m)g2κaη, (27) ξ=11-β-mBg2κa1-(3mB-m)g2κaη, (27)

v=ξω2π. (28) v=ξω2π. (28)

As the duty factor β increases, ξ and v increase.

When k = 2 and ψˆ1≈2βπ ψ^1≈2βπ , the swing movements of the solution propagate from anterior to posterior with the wavelength from (27) and velocity from (28). This is similar to the retrograde-wave gait of centipedes (Full, 1997; Kuroda et al., 2014). In contrast, when ψˆN≈2(1−β)π ψ^N≈2(1-β)π in (24), the solution, whose front and rear parts are reversed from the solution of the retrograde-wave gait, is similar to the direct-wave gait of millipedes (Full, 1997; Kuroda et al., 2014). Our model has solutions corresponding to direct-, retrograde-, and source-wave gaits for N ≥ 3.

4. Discussion

4.1. Body Elasticity and Local Sensory Feedback Generate Interlimb Coordination

We assume that animal gaits can be represented as attractors of dynamic systems, as in Schöner et al. (1990), and that the essential structure of the dynamic system can be extracted by using a simple model (Full and Koditschek, 1999). Accumulating an understanding of such simple models allows us to understand actual complex phenomena. In the present study, we focused on the embodied sensorimotor interaction to generate multi-legged locomotion in a decentralized manner. We showed that the local sensory feedback and phase resetting generates the direct-, retrograde-, and source-wave gaits observed in multi-legged animals using a simple model to extract the essential features. We found that body elasticity is a key to generating the interlimb coordination. Specifically, the body natural frequency must be larger than the gait frequency ( κ/mi−−−−√≫ω κ/mi≫ω ). The analytical representation showed the existence of these gaits for N ≥ 3, and revealed the parameter domain as follows: the boundary condition of the mass, m/2 < m_B < m, and the amplitude of the swing leg movement, mg/κ < a < L. These findings improve our understanding of the mechanism of interlimb coordination in multi-legged locomotion.

So far, many studies have investigated the effects of sensory feedback on interlimb coordination. Owaki et al. (2012), Fukuoka et al. (2015), and Owaki and Ishiguro (2017) showed that quadruped robots, whose legs are controlled by distributed oscillators with load sensory feedback, generate walking, trotting, and galloping gaits depending on the speed. Tamura et al. (2016), Kano et al. (2017), and Yasui et al. (2017) also showed that load sensory feedback generates a millipede-like direct-wave gait. Our previous studies showed that phase resetting induces gait transitions in quadruped locomotion (Aoi et al., 2011, 2013) and hexapod locomotion (Fujiki et al., 2013; Ambe et al., 2015, 2018). In the present study, we also demonstrated that phase resetting generates retrograde-, source-, and direct-wave gaits in multi-legged locomotion. These findings are due to the embodied sensorimotor interaction, which is critical for the interlimb coordination regardless of the number of legs. Our analytical description using a simple model is helpful for clarifying the functional roles of the embodied sensorimotor interaction in interlimb coordination.

4.2. Relation With Intersegmental Coordination in Other Organisms

Other organisms also show intersegmental coordination in their locomotion, which appears as waves. For example, insect larvae and worms exhibit direct waves while crawling (Trimmer and Issberner, 2007; Paoletti and Mahadevan, 2014). Lampreys, leeches, and roundworms (C. Elegans) show body undulations that propagate backward during swimming. Motile cilia show metachronal waves. We discuss the relationship between our results and such intersegmental coordination by focusing on local sensory feedback.

Paoletti and Mahadevan (2014) developed an earthworm crawling model using mass points connected by spring dampers in an asymmetric friction environment. The model is controlled by neuro-muscular dynamics with the local sensory feedback, which contracts the segment muscles based on a stretch threshold and generates metachronal waves. Umedachi et al. (2016) proposed a similar model to reproduce the direct wave crawling motion of larvae in a decentralized manner. Each segment repeats stretching and contraction with distributed CPGs with local sensory feedback of the friction and velocity information for the segment. These results are similar to ours in the sense that metachronal waves are generated by local sensory feedback. However, their coordination is generated by asymmetric environmental friction, whereas ours is by different foot-contact timings.

The undulation patterns of lampreys and leeches as they swim are mainly generated by internal coupling of CPGs (Cohen et al., 1992; Grillner et al., 1995; KristanJr. et al., 2005). However, intersegmental coordination in leeches is not disrupted much even if the ganglion is cut (Yu et al., 1999), which suggests the importance of sensory feedback for segmental coordination. In addition, there is no clear evidence for the presence of CPGs in C. Elegans, and it has been suggested that proprioceptive receptors play a significant role in coordination (Wen et al., 2012). Boyle et al. (2012) modeled C. Elegans as a series of simple segments connected by elastic elements. They reproduced body undulations by neuro-muscular dynamics without CPGs. The dynamics of each segment is affected by the sensory feedback, which integrates several posterior segment stretches. When the feedback has only local interaction, such as stretch of the segment, no coordinated wave appears in a less viscous environment, such as water, and it appears only in a highly viscous fluid environment. That is, the fluid viscosity is responsible for generating the coordinated wave.

Motile cilia of organisms generate metachronal waves in a decentralized manner due to local interaction between the environmental fluid and the cilia (Elgeti and Gompper, 2013). Each cilium is controlled by simple switching inputs of power stroke and return stroke. These inputs switch when the cilium achieves a certain curvature; that is, it receives local sensory feedback about the curvature. Although each cilium moves independently, the metachronal wave for cilia motion appears through local interaction by the flow of the neighboring fluid.

An analytical description would help to understand these intersegmental coordination mechanisms. Thus, in the future, we would like to develop simple physical models for intersegmental coordination.

4.3. Physical Explanation of Countless Solutions

Our model has countless solutions as obtained by the serially connected set of the retrograde-, source-, and direct-wave gaits (Figure 11). This non-uniqueness is because our model has a conservative quantity as explained below.

We define E^j on the Poincaré section Σ^j (j ∈ [1, N + 1]) by

Ej=∑i=1Nψji. (29) Ej=∑i=1Nψij. (29)

From (11), the Poincaré map of E^j is represented as

Ej↦Ej+ϕ11−ϕN+1N+1. (30) Ej↦Ej+ϕ11-ϕN+1N+1. (30)

The relation ϕ11=ϕN+1N+1=ϕˆtd ϕ11=ϕN+1N+1=ϕ^td holds if the neighboring legs (Legs 2 and N) are in the stance phases ( ϕ12∈[0,2βπ] ϕ21∈[0,2βπ] and ϕN+1N∈[0,2βπ] ϕNN+1∈[0,2βπ] ) when each leg (Legs 1 and N + 1) touches the ground (this relation is also satisfied for the solution sets). In this case, (30) becomes E^j ↦ E^j, where E^j is a conservative quantity, which produces countless solutions.

To explain the physical meaning of this conservative quantity, we represent ϕ11 ϕ11 and ϕN+1N+1 ϕN+1N+1 using (3) and (12) when ϕ12∈[0,2βπ] ϕ21∈[0,2βπ] and ϕN+1N∈[0,2βπ] ϕNN+1∈[0,2βπ] .

ϕ11=2π−(1−β)m1gπκa, (31) ϕ11=2π-(1-β)m1gπκa, (31)

ϕN+1N+1=2π−(1−β)mN+1gπκa. (32) ϕN+1N+1=2π-(1-β)mN+1gπκa. (32)

These equations show that ϕ11=ϕN+1N+1 ϕ11=ϕN+1N+1 is satisfied only when m₁ = m_{N + 1}. Thus, this conservative quantity is derived from the symmetry between the front and tail masses (m₁ = m_{N + 1} = m_B). Strictly speaking, the symmetry of other parameters, such as κ and a, is also required for the conservative quantity.

Although our model generates multiple gaits because of the symmetry, animals prefer specific gaits. For example, centipedes prefer the retrograde-wave gait (Full, 1997; Kuroda et al., 2014). In contrast, millipedes prefer the direct-wave gait to move forward (Full, 1997; Kuroda et al., 2014) and use the source-wave gait when the body axis is bent like a U shape (Tamura et al., 2016). Our model does not explain these preferences. However, when the symmetry is broken, our model has no conservative quantity and generates a specific gait depending on the parameters. We would like to investigate the gait preference by incorporating asymmetries into our model in future studies.

4.4. Limitation of This Study

This study demonstrated that phase resetting generates direct-, retrograde-, and source-wave gaits in an elastic body whose natural frequency is larger than the gait frequency. However, our model did not replicate the details of animal gaits. For example, our model did not explain the gait preference of animals, as described above. While centipedes and millipedes use different contralateral interlimb coordination (i.e., left–right antiphase and in-phase movements, respectively), our model did not consider this behavior. Furthermore, although the swinging of multiple legs propagates simultaneously in millipedes and centipedes (Kuroda et al., 2014), only one leg swing is propagated in our model. These limitations and discrepancies come from the simple way in which our model extracts the essence of multi-leg motion, which was the focus in the present study. However, our results give clues to overcome these limitations and discrepancies. For example, the gait preference may be explained by introducing an asymmetry into our model, such as different masses between the head and tail, to eliminate the conservative quantity in our model. Better understanding of a simple model will provide insights for the design of a more complicated model. The extension of our model is one of our essential future tasks.

4.5. Application to Legged Robots

Multi-legged robots have been developed to extract the essence of dexterous traveling ability. On one hand, the flexible body axis contributes to rapid movement or ability to traverse rough terrain. In previous work (Aoi et al., 2016), we showed that high body-axis flexibility induces body undulations through Hopf bifurcation, which contributes to rapid turning. Hoffman and Wood (2011) showed that a passive undulatory gait increases the locomotion speed. Koh et al. (2010) and Kinugasa et al. (2017) demonstrated that multi-legged robots with a flexible body axis show high mobility in uneven terrain.

On the other hand, contacting the ground with many legs is useful for avoiding stumbling in various environments. Inagaki et al. (2010) proposed a distributed control method in which the legs follow the contact points of anterior legs, which allowed a robot to walk in various environments as long as the front legs choose solid footholds. Hayakawa et al. (2020) proposed a gait generation strategy to ensure static stability for single-legged modular robots to create a cluster with various leg configurations. Kano et al. (2017) showed the adaptability when part of the terrain is removed.

For static gaits, sensory feedback is useful for generating interlimb coordination, as in the present study, and to gain better adaptability to the environment (Kano et al., 2017; Yasui et al., 2017). To design a method for controlling multi-legged robots with sensory feedback, it is important to understand how and when the sensory feedback affects walking motion. The present study focused on phase resetting at foot contact as the sensory feedback and clarified the effects by using the analytical description of a simple model. We showed that phase resetting contributes to generating the coordinated gaits in a decentralized manner via embodied sensorimotor interaction when the model has an elastic body axis whose natural frequency is larger than the gait frequency. We also analytically derived the range of physical parameters where our analysis is valid. This helps when designing a controller to generate various gaits for multi-legged robots. In particular, when we design the leg movement using oscillator phases, as in Aoi et al. (2017) and Ambe et al. (2018), phase resetting will generate coordinated gaits regardless of the number of legs. We would like to investigate it through robot experiments in future studies.

5. Conclusion

This study used a simple model to analytically reveal that local sensory feedback, phase resetting, generates the direct-, retrograde-, and source-wave gaits observed in multi-legged animals in a decentralized manner via the embodied sensorimotor interaction. The model comprises massless legs and mass points connected by vertical springs to imitate a flexible body. The phase oscillators control the vertical movements of the legs independently, and the phases are only affected by phase resetting upon foot contact. The dynamic simulations show that countless periodic solutions of the three gaits emerge depending on the initial phase. Furthermore, these analytical solutions were derived under some assumptions deduced from the simulation results. They showed that all three gaits exist regardless the number of legs and revealed the parameter domain. The reason for the coexistence of the three gaits is explained by a conservative quantity due to the front–rear symmetry of our model. Because our model is limited to specific situations, such as the front–rear symmetry, permitting only vertical movements while ignoring the contralateral interlimb coordination, and propagation of only one swing leg within one wavelength, we need to incorporate more realistic situations into our model in the future.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

YA conceived of the study, designed the study, carried out the simulation and theoretical analysis, and drafted the manuscript. SA helped conceiving and designing of the study, improved the theoretical analysis, and critically revised the manuscript. KT helped conceiving and designing of the study and critically revised the manuscript. FM helped conceiving and designing of the study, coordinated the study, and critically revised the manuscript. All authors gave final approval for publication and agree to be held accountable for the work performed therein.

Funding

This work was supported by JSPS KAKENHI Grant Number JP18K13723.

Conflict of Interest

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.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References

Adachi, M., Aoi, S., Kamimura, T., Tsuchiya, K., and Matsuno, F. (2020). Body torsional flexibility effects on stability during trotting and pacing based on a simple analytical model. Bioinspir. Biomim. 15:055001. doi: 10.1088/1748-3190/ab968d

Akay, T., Haehn, S., Schmitz, J., and Büschges, A. (2004). Signals from load sensors underlie interjoint coordination during stepping movements of the stick insect leg. J. Neurophysiol. 92, 42–51. doi: 10.1152/jn.01271.2003

Ambe, Y., Aoi, S., Nachstedt, T., Manoonpong, P., Wörgötter, F., and Matsuno, F. (2018). Simple analytical model reveals the functional role of embodied sensorimotor interaction in hexapod gaits. PLoS ONE, Kyoto, 13:e0192469. doi: 10.1371/journal.pone.0192469

Ambe, Y., Nachstedt, T., Manoonpong, P., Wörgötter, F., and Matsuno, F. (2015). “Embodied sensorimotor interaction for insect locomotion,” in The First International Symposium on Swarm Behavior and Bio-Inspired Robotics, 340–343.

Aoi, S., Katayama, D., Fujiki, S., Tomita, N., Funato, T., Yamashita, T., et al. (2013). A stability-based mechanism for hysteresis in the walk-trot transition in quadruped locomotion. J. R. Soc. Interface 10:20120908. doi: 10.1098/rsif.2012.0908

Aoi, S., Manoonpong, P., Ambe, Y., Matsuno, F., and Wörgötter, F. (2017). Adaptive control strategies for interlimb coordination in legged robots: a review. Front. Neurorobot. 11:39. doi: 10.3389/fnbot.2017.00039

Aoi, S., Tanaka, T., Fujiki, S., Funato, T., Senda, K., and Tsuchiya, K. (2016). Advantage of straight walk instability in turning maneuver of multilegged locomotion: a robotics approach. Sci. Rep. 6:30199. doi: 10.1038/srep30199

Aoi, S., Yamashita, T., and Tsuchiya, K. (2011). Hysteresis in the gait transition of a quadruped investigated using simple body mechanical and oscillator network models. Phys. Rev. E 83:061909. doi: 10.1103/PhysRevE.83.061909

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

Büschges, A. (1995). Role of local nonspiking interneurons in the generation of rhythmic motor activity in the stick insect. J. Neurobiol. 27, 488–512. doi: 10.1002/neu.480270405

Büschges, A., Akay, T., Gabriel, J. P., and Schmidt, J. (2008). Organizing network action for locomotion: Insights from studying insect walking. Brain Res. Rev. 57, 162–171. doi: 10.1016/j.brainresrev.2007.06.028

Cohen, A. H., Ermentrout, G. B., Kiemel, T., Kopell, N., Sigvardt, K. A., and Williams, T. L. (1992). Modelling of intersegmental coordination in the lamprey central pattern generator for locomotion. Trends Neurosci. 15, 434–438. doi: 10.1016/0166-2236(92)90006-T

Cruse, H., Kindermann, T., Schumm, M., Dean, J., and Schmitz, J. (1998). Walknet-a biologically inspired network to control six-legged walking. Neural Netw. 11, 1435–1447. doi: 10.1016/S0893-6080(98)00067-7

Daun, S., Rubin, J., and Rybak, I. (2009). Control of oscillation periods and phase durations in half-center central pattern generators: a comparative mechanistic analysis. J. Comput. Neurosci. 27, 3–36. doi: 10.1007/s10827-008-0124-4

Daun-Gruhn, S., and Büschges, A. (2011). From neuron to behavior: dynamic equation-based prediction of biological processes in motor control. Biol. Cybern. 105, 71–88. doi: 10.1007/s00422-011-0446-6

Delcomyn, F. (2004). Insect walking and robotics. Annu. Rev. Entomol. 49, 51–70. doi: 10.1146/annurev.ento.49.061802.123257

Elgeti, J., and Gompper, G. (2013). Emergence of metachronal waves in cilia arrays. Proc. Natl. Acad. Sci. U.S.A. 110, 4470–4475. doi: 10.1073/pnas.1218869110

Fujiki, S., Aoi, S., Funato, T., Tomita, N., Senda, K., and Tsuchiya, K. (2013). Hysteresis in the metachronal-tripod gait transition of insects: a modeling study. Phys. Rev. E 88:012717. doi: 10.1103/PhysRevE.88.012717

Fukuoka, Y., Habu, Y., and Fukui, T. (2015). A simple rule for quadrupedal gait generation determined by leg loading feedback: a modeling study. Sci. Rep. 5:8169. doi: 10.1038/srep08169

Full, R. J. (1997). “Invertebrate locomotor system,” in The Handbook of Comparative Physiology, ed W. Dantzler (Oxford: Oxford University Press), 853–930. doi: 10.1002/cphy.cp130212

Full, R. J., and Koditschek, D. E. (1999). Templates and anchors: neuromechanical hypotheses of legged locomotion on land. J. Exp. Biol. 202, 3325–3332. doi: 10.1242/jeb.202.23.3325

Grillner, S., Deliagina, T., el Manira, A., Hill, R. H., Lansner, A., Orlovsky, G. N., et al. (1995). Neural networks that co-ordinate locomotion and body orientation in lamprey. Trends Neurosci. 18, 270–279. doi: 10.1016/0166-2236(95)80008-P

Hayakawa, T., Kamimura, T., Kaji, S., and Matsuno, F. (2020). Autonomous distributed system for gait generation for single-legged modular robots connected in various configurations. IEEE Trans. Robot. 36, 1491–1510. doi: 10.1109/TRO.2020.2992983

Hoffman, K. L., and Wood, R. J. (2011). “Passive undulatory gaits enhance walking in a myriapod millirobot,” in 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems (San Francisco, CA), 1479–1486. doi: 10.1109/IROS.2011.6048268

Hooper, S. L., Guschlbauer, C., Blümel, M., Rosenbaum, P., Gruhn, M., Akay, T., et al. (2009). Neural control of unloaded leg posture and of leg swing in stick insect, cockroach, and mouse differs from that in larger animals. J. Neurosci. 29, 4109–4119. doi: 10.1523/JNEUROSCI.5510-08.2009

Ijspeert, A. J. (2001). A connectionist central pattern generator for theaquatic and terrestrial gaits of a simulated salamander. Biol. Cybern. 84, 331–348. doi: 10.1007/s004220000211

Ijspeert, A. J. (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

Inagaki, S., Niwa, T., and Suzuki, T. (2010). “Follow-the-contact-point gait control of centipede-like multi-legged robot to navigate and walk on uneven terrain,” in 2010 IEEE/RSJ International Conference on Intelligent Robots and Systems (Taipei), 5341–5346. doi: 10.1109/IROS.2010.5651324

Kamimura, T., Aoi, S., Higurashi, Y., Wada, N., Tsuchiya, K., and Matsuno, F. (2021). Dynamical determinants enabling two different types of flight in cheetah gallop to enhance speed through spine movement. Sci. Rep. 11:9631. doi: 10.1038/s41598-021-88879-0

Kano, T., Sakai, K., Yasui, K., Owaki, D., and Ishiguro, A. (2017). Decentralized control mechanism underlying interlimb coordination of millipedes. Bioinspir. Biomim. 12:036007. doi: 10.1088/1748-3190/aa64a5

Kinugasa, T., Osuka, K., Hayashi, R., Miyamoto, N., and Yoshida, K. (2017). Development of a small and lightweight myriapod robot using passive dynamics. Artif. Life Robot. 22, 429–434. doi: 10.1007/s10015-017-0378-x

Koh, D., Yang, J., and Kim, S. (2010). “Centipede robot for uneven terrain exploration: design and experiment of the flexible biomimetic robot mechanism,” in 2010 3rd IEEE RAS EMBS International Conference on Biomedical Robotics and Biomechatronics (Tokyo), 877–881. doi: 10.1109/BIOROB.2010.5627776

Kristan, W. B. Jr, Calabrese, R. L., and Friesen, W. O. (2005). Neuronal control of leech behavior. Prog. Neurobiol. 76, 279–327. doi: 10.1016/j.pneurobio.2005.09.004

Kuroda, S., Kunita, I., Tanaka, Y., Ishiguro, A., Kobayashi, R., and Nakagaki, T. (2014). Common mechanics of mode switching in locomotion of limbless and legged animals. J. R. Soc. Interface 11:20140205. doi: 10.1098/rsif.2014.0205

Matsuoka, K. (1987). Mechanisms of frequency and pattern control inthe neural rhythm generators. Biol. Cybern. 56, 345–353. doi: 10.1007/BF00319514

Owaki, D., and Ishiguro, A. (2017). A quadruped robot exhibiting spontaneous gait transitions from walking to trotting to galloping. Sci. Rep. 7:277. doi: 10.1038/s41598-017-00348-9

Owaki, D., Kano, T., Nagasawa, K., Tero, A., and Ishiguro, A. (2012). Simple robot suggests physical interlimb communication is essential for quadruped walking. J. R. Soc. Interface 10:20120669. doi: 10.1098/rsif.2012.0669

Paoletti, P., and Mahadevan, L. (2014). A proprioceptive neuromechanical theory of crawling. Proc. R. Soc. B 281:20141092. doi: 10.1098/rspb.2014.1092

Pasemann, F., Hild, M., and Zahedi, K. (2003). “SO(2)-networks as neural oscillators,” in IWANN 2003: Computational Methods in Neural Modeling, Vol. 2686, eds J. Mira and J. R. Álvarez (Springer), 141–151. doi: 10.1007/3-540-44868-3_19

Proctor, J. L., and Holmes, P. (2018). The effects of feedback on stability and maneuverability of aphase-reduced model for cockroach locomotion. Biol. Cybern. 112, 387–401. doi: 10.1007/s00422-018-0762-1

Rybak, I. A., Shevtsova, N. A., Lafreniere-Roula, M., and McCrea, D. A. (2006). Modelling spinal circuitry involved in locomotor patterngeneration: insights from deletions during fictivelocomotion. J. Physiol. 577, 617–639. doi: 10.1113/jphysiol.2006.118703

Schilling, M., Hoinville, T., Schmitz, J., and Cruse, H. (2013). Walknet, a bio-inspired controller for hexapod walking. Biol. Cybern. 107, 397–419. doi: 10.1007/s00422-013-0563-5

Schomburg, E. D., Petersen, N., Barajon, I., and Hultborn, H. (1998). Flexor reflex afferents reset the step cycle during fictive locomotion in the cat. Exp. Brain Res. 122, 339–350. doi: 10.1007/s002210050522

Schöner, G., Jiang, W., and Kelso, J. (1990). A synergetic theory of quadrupedal gaits and gait transitions. J. Theor. Biol. 142, 359–391. doi: 10.1016/S0022-5193(05)80558-2

Seipel, J. E., Holmes, P. J., and Full, R. J. (2004). Dynamics and stability of insect locomotion: a hexapedal model for horizontal plane motions. Biol. Cybern. 91, 76–90. doi: 10.1007/s00422-004-0498-y

Shaw, K. M., Lyttle, D. N., Gill, J. P., Cullins, M. J., McManus, J. M., Lu, H., et al. (2015). The significance of dynamical architecture for adaptive responses to mechanical loads during rhythmic behavior. J. Comput. Neurosci. 38, 25–51. doi: 10.1007/s10827-014-0519-3

Spence, A. J., Revzen, S., Seipel, J., Mullens, C., and Full, R. J. (2010). Insects running on elastic surfaces. J. Exp. Biol. 213, 1907–1920. doi: 10.1242/jeb.042515

Steingrube, S., Timme, M., Wörgötter, F., and Manoonpong, P. (2010). Self-organized adaptation of a simple neural circuit enables complex robot behaviour. Nat. Phys. 6, 224–230. doi: 10.1038/nphys1508

Tamura, T., Satou, E., and Inoue, K. (2016). “Constructive understanding of decentralized leg motion control mechanism of myriapods,” in 28th Sympsium of Destributed Autonomous System, Hiroshima, 149–152. (in Japanese).

Tanase, M., Ambe, Y., Aoi, S., and Matsuno, F. (2015). A galloping quadruped model using left-right asymmetry in touchdown angles. J. Biomech. 48, 3383–3389. doi: 10.1016/j.jbiomech.2015.06.003

Ting, L. H., and Chiel, H. J. (2017). “Chapter 12: Muscle, biomechanics, and implications for neural control,” in Neurobiology of Motor Control: Fundamental Concepts and New Directions, eds S. L. Hooper and A. Büschges (Hoboken, NJ: Wiley), 365–416. doi: 10.1002/9781118873397.ch12

Trimmer, B., and Issberner, J. (2007). Kinematics of soft-bodied, legged locomotion in manduca sexta larvae. Biol. Bull. 212, 130–142. doi: 10.2307/25066590

Tuthill, J. C., and Wilson, R. I. (2016). Mechanosensation and adaptive motor control in insects. Curr. Biol. 26, 1022–1038. doi: 10.1016/j.cub.2016.06.070

Umedachi, T., Kano, T., Ishiguro, A., and Trimmer, B. A. (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

van der Pol, B. (1926). On relaxation oscillations. Lond. Edinb. Dublin Philos. Mag. J. Sci. 2, 978–992. doi: 10.1080/14786442608564127

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

Yasui, K., Kano, T., Standen, E. M., Aonuma, H., Ijspeert, A. J., and Ishiguro, A. (2019). Decoding the essential interplay between central and peripheral control in adaptive locomotion of amphibious centipedes. Sci. Rep. 9:18288. doi: 10.1038/s41598-019-53258-3

Yasui, K., Sakai, K., Kano, T., Owaki, D., and Ishiguro, A. (2017). Decentralized control scheme for myriapod robot inspired by adaptive and resilient centipede locomotion. PLoS ONE 12:e171421. doi: 10.1371/journal.pone.0171421

Yu, X., Nguyen, B., and Friesen, W. O. (1999). Sensory feedback can coordinate the swimming activity of the leech. J. Neurosci. 19, 4634–4643. doi: 10.1523/JNEUROSCI.19-11-04634.1999