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

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

$\begin{array}{l}{m}_i{\ddot{x}}_i=\\ \{\begin{array}{cc}\begin{array}{l}-\kappa (x_1-x_2)-\sigma ({\dot{x}}_1-{\dot{x}}_2)-m_{1}g\hfill \\ -p_1\left\{K(x_1-l_1)+D({\dot{x}}_1-{\dot{l}}_1)\right\}\hfill \end{array}& i=1\\ \begin{array}{l}-\kappa (2x_i-x_{i+1}-x_{i-1})-\sigma (2{\dot{x}}_i-{\dot{x}}_{i+1}-{\dot{x}}_{i-1})-m_{i}g\hfill \\ -p_i\left\{K(x_i-l)+D({\dot{x}}_i-{\dot{l}}_i)\right\}\hfill \end{array}& i=2,\dots ,N\\ \begin{array}{l}-\kappa (x_{N+1}-x_N)-\sigma ({\dot{x}}_{N+1}-{\dot{x}}_N)-m_{N+1}g\hfill \\ -p_{N+1}\left\{K(x_{N+1}-l_{N+1})+D({\dot{x}}_{N+1}-{\dot{l}}_{N+1})\right\}\hfill \end{array}& i=N+1\end{array}& (1)\end{array}$

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

$\begin{array}{ll}{p}_i=\{\begin{array}{cc}0& l_i<x_i\\ 1& \text{otherwise}.\end{array}& (2)\end{array}$

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

$\begin{array}{ll}{l}_i=l({\varphi }_i)=\{\begin{array}{cc}L& 0\le {\varphi }_i\le 2\beta \pi \\ L-a\frac{{\varphi }_i-2\beta \pi }{(1-\beta )\pi }& 2\beta \pi <{\varphi }_i<(1+\beta )\pi \\ L-a\frac{2\pi -{\varphi }_i}{(1-\beta )\pi }& (1+\beta )\pi \le {\varphi }_i<2\pi \end{array}& (3)\end{array}$

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:

$\begin{array}{ll}\frac{d{\varphi }_i}{dt}=\omega +y_i,& (4)\end{array}$

where ω is the basic frequency of locomotion and y_i represents a sensory feedback. We used $\omega \ll \sqrt{\kappa /m_i}$ to make the model walk slowly.

FIGURE 2

Figure 2. (A) Phase oscillator model. (B) Leg length l_i depending on phase ϕ_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

$\begin{array}{ll}\tau \frac{d{y}_i}{dt}=-y_i+u_i& (5)\end{array}$

$\begin{array}{ll}{u}_i=\{\begin{array}{cc}0& 0\le {\varphi }_i<(1+\beta )\pi \\ (2\pi -{\varphi }_i)\delta (t-t_i^i)& (1+\beta )\pi \le {\varphi }_i<2\pi \end{array}& (6)\end{array}$

where $t_i^i$ 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

$\begin{array}{ll}{\psi }_i={\varphi }_{i+1}-{\varphi }_i.& (7)\end{array}$

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 ${\Psi }^1\equiv {\left[{\psi }_1^1{\psi }_2^1{\psi }_3^1\right]}^{\text{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 Ψ¹ simulated from six sets of initial values: (A) ${\psi }_1^1$, (B) ${\psi }_2^1$, and (C) ${\psi }_3^1$.

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 ${\psi }_1^1$–${\psi }_2^1$ and ${\psi }_2^1$–${\psi }_3^1$ 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 Ψ¹ of periodic solutions: (A) shown in ${\psi }_1^1$-${\psi }_2^1$-${\psi }_3^1$ space, (B) projected to a plane parallel to the ${\psi }_1^1$-${\psi }_2^1$ plane, and (C) projected to a plane parallel to the ${\psi }_2^1$-${\psi }_3^1$ 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 ${\psi }_3^1=2(1-\beta )\pi ,\pi$, 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) ${\psi }_3^1=2(1-\beta )\pi$, (B) ${\psi }_3^1=\pi$, and (C) ${\psi }_3^1=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:

$\begin{array}{l}{\widehat{\psi }}_1^1\approx 2(1-\beta )\pi ,{\widehat{\psi }}_2^1\approx 2(1-\beta )\pi ,2(1-\beta )\pi ⪅{\widehat{\psi }}_3^1⪅2\beta \pi \end{array}$

where $\widehat{()}$ indicates the periodic solution. This gait appears as one line segment, as shown in Figure 4. Figure 6A shows the footprint diagram for ${\widehat{\Psi }}^1={\left[1.130.981.50\right]}^{\text{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:

$\begin{array}{l}2(1-\beta )\pi ⪅{\widehat{\psi }}_1^1⪅2\beta \pi ,{\widehat{\psi }}_2^1\approx 2\beta \pi ,{\widehat{\psi }}_3^1\approx 2\beta \pi .\end{array}$

This also appears as one line segment (to be precise, it is slightly curved) in Figure 4. Figure 6B shows the footprint diagram for ${\widehat{\Psi }}^1={\left[3.975.295.29\right]}^{\text{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:

$\begin{array}{l}{\widehat{\psi }}_1^1\approx 2(1-\beta )\pi ,2(1-\beta )\pi ⪅{\widehat{\psi }}_2^1⪅2\beta \pi ,{\widehat{\psi }}_3^1\approx 2\beta \pi .\end{array}$

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 ∈ [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 ${\varphi }_i^i$. These assumptions give a touchdown relation between ${\varphi }_i^i,{\varphi }_{i-1}^i$, and ${\varphi }_{i+1}^i$. 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 ${\varphi }_{i-1}^{i-1},{\varphi }_i^i,{\varphi }_{i+1}^{i+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 $\omega \ll \sqrt{\kappa /m_i}$ 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

$\begin{array}{ll}\{\begin{array}{cc}\kappa (x_2-x_1)+R_1=m_{1}g& i=1\\ \kappa (x_{i-1}-x_i)+\kappa (x_{i+1}-x_i)+R_i=m_{i}g& i\in \left[2,N\right]\\ \kappa (x_N-x_{N+1})+R_{N+1}=m_{N+1}g& i=N+1\end{array},& (8)\end{array}$

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

$\begin{array}{ll}\frac{d{\varphi }_i}{dt}=\omega +u_i.& (9)\end{array}$

This gives

$\begin{array}{ll}\frac{d{\psi }_i}{dt}=u_{i+1}-u_{i}i=1,\dots ,N.& (10)\end{array}$

A Poincaré section Σ^j (j = 1, …, N+1) was defined for the relative phases just before the touchdown of Leg j and we used ${\Psi }^j\equiv {\left[{\psi }_1^j{\psi }_2^j\dots {\psi }_N^j\right]}^{\text{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

$\begin{array}{ll}{\psi }_i^j\mapsto {\psi }_i^j+{\varphi }_i^i-{\varphi }_{i+1}^{i+1}i=1,\dots ,N.& (11)\end{array}$

To obtain the Poincaré map, ${\varphi }_i^i(i=1,\dots ,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, $x_{i-1}=l({\varphi }_{i-1}^i)$, $x_i=l({\varphi }_i^i)$, $x_{i+1}=l({\varphi }_{i+1}^i)$, 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

$\begin{array}{ll}l({\varphi }_i^i)=\{\begin{array}{cc}{f}_{\text{B}}({\varphi }_2^i)& i=1\\ f({\varphi }_{i-1}^i,{\varphi }_{i+1}^i)& i\in \left[2,N\right]\\ f_{\text{B}}({\varphi }_N^i)& i=N+1\\ \end{array},& (12)\end{array}$

where

$\begin{array}{ll}f({\varphi }_{\text{h}},{\varphi }_{\text{f}})=\frac{1}{2}\left\{l({\varphi }_{\text{h}})+l({\varphi }_{\text{f}})\right\}-\frac{mg}{2\kappa }& (13)\end{array}$

$\begin{array}{l}{f}_{\text{B}}(\varphi )=l(\varphi )-\frac{m_{\text{B}}g}{\kappa }.\end{array}$

The substitution of (12) into (3) gives one relation between ${\varphi }_{i-1}^i,{\varphi }_i^i$, and ${\varphi }_{i+1}^i$ for each i (i = 1, …, N + 1). Because we can write ${\varphi }_{i-1}^i$ and ${\varphi }_{i+1}^i$ as ${\varphi }_{i-1}^{i-1},{\varphi }_i^i,{\varphi }_{i+1}^{i+1}$, and Ψ^j, as discussed in the following section, ${\varphi }_i^i$ 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 ${\varphi }_i^i$ (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 ${\varphi }_i^i$ by accounting for the effects of phase resetting. Second, we derive the conditions of ${\varphi }_i^i$ 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), ${\varphi }_{k-1}^k$ is in the stance or takeoff phase ($0\le {\varphi }_{k-1}^k<(1+\beta )\pi$), and ${\varphi }_{k+1}^k$ is in the takeoff phase ($2\beta \pi \le {\varphi }_{k+1}^k<(1+\beta )\pi$). For the touchdown of Leg j ∈ [2, k − 1] (direct-wave region), ${\varphi }_{j-1}^j$ is in the takeoff phase ($2\beta \pi \le {\varphi }_{j-1}^j<(1+\beta )\pi$) and ${\varphi }_{j+1}^j$ is in the stance phase ($0\le {\varphi }_{j+1}^j<2\beta \pi$) because the swing movements of Legs 1 to k − 1 propagate anteriorly. For the touchdown of Leg i ∈ [k + 1, N] (retrograde-wave region), ${\varphi }_{i-1}^i$ is in the stance phase ($0\le {\varphi }_{i-1}^i<2\beta \pi$) and ${\varphi }_{i+1}^i$ is in the takeoff phase ($2\beta \pi \le {\varphi }_{i+1}^i<(1+\beta )\pi$) 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), ${\varphi }_2^1$ and ${\varphi }_N^{N+1}$ are in the stance phases ($0\le {\varphi }_2^1,{\varphi }_N^{N+1}<2\beta \pi$).

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 ${\varphi }_n^k\in ({\varphi }_n^{n-1},{\varphi }_n^n]$, 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. ${\psi }_{i-1}^i$ and ${\psi }_i^i$ just before the touchdown of Leg i (i = k, n, k+1, k+2, …, n−2, n−1) are represented by Ψ^k, ${\varphi }_{i-1}^{i-1}$, and ${\varphi }_{i+1}^{i+1}$, as highlighted in orange.

When the order of touchdown events is determined, the relative phases ${\psi }_{i-1}^i$ and ${\psi }_i^i$ just before touchdown of Leg i can be represented using Ψ^k and ${\varphi }_{i-1}^{i-1}$, and ${\varphi }_{i+1}^{i+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, ${\psi }_{i-1}^i$ and ${\psi }_i^i$ just before touchdown of Leg i can be represented using Ψ^k, ${\varphi }_{i-1}^{i-1}$, and ${\varphi }_{i+1}^{i+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, ${\psi }_{i-1}^i$ and ${\psi }_i^i$ are given by ${\psi }_{i-1}^i={\psi }_{i-1}^k$ and ${\psi }_i^i={\psi }_i^k$, 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, ${\varphi }_{i-1}^i={\varphi }_i^i-{\psi }_{i-1}^k$ and ${\varphi }_{i+1}^i={\varphi }_i^i+{\psi }_i^k-2\pi$ (${\varphi }_{i+1}^i$ is subtracted by 2π to satisfy $0\le {\varphi }_{i+1}^i<2\pi$) are satisfied because ${\psi }_{i-1}^i={\varphi }_i^i-{\varphi }_{i-1}^i$ and ${\psi }_i^i={\varphi }_{i+1}^i-{\varphi }_i^i$. For i ∈ [k+1, n−2], ${\psi }_{i-1}^i$ and ${\psi }_i^i$ are given by ${\psi }_{i-1}^i={\psi }_{i-1}^k-(2\pi -{\varphi }_{i-1}^{i-1})$ and ${\psi }_i^i={\psi }_i^k$, 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, ${\varphi }_{i-1}^i={\varphi }_i^i-{\psi }_{i-1}^k+2\pi -{\varphi }_{i-1}^{i-1}$ and ${\varphi }_{i+1}^i={\varphi }_i^i+{\psi }_i^k-2\pi$ are satisfied. For i = n−1, ${\psi }_{i-1}^i$ and ${\psi }_i^i$ are given by ${\psi }_{i-1}^i={\psi }_{i-1}^k-(2\pi -{\varphi }_{i-1}^{i-1})$ and ${\psi }_i^i={\psi }_i^k+(2\pi -{\varphi }_{i+1}^{i+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, ${\varphi }_{i-1}^i={\varphi }_i^i-{\psi }_{i-1}^k+2\pi -{\varphi }_{i-1}^{i-1}$ and ${\varphi }_{i+1}^i={\varphi }_i^i+{\psi }_i^k-{\varphi }_{i+1}^{i+1}$ are satisfied. This means that ${\varphi }_{i-1}^i$ and ${\varphi }_{i+1}^i$ can be represented by Ψ^k, ${\varphi }_{i-1}^{i-1}$, ${\varphi }_i^i$, and ${\varphi }_{i+1}^{i+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 ${\psi }_i^k(i=1,\dots ,N)$ and ${\varphi }_i^i(i=1,\dots ,N+1)$ as

$\begin{array}{ll}l({\varphi }_i^i)=\{\begin{array}{cc}{f}_{\text{B}}({\varphi }_1^1+{\psi }_1^k-2\pi )& i\in S_1^k\\ f_{\text{B}}({\varphi }_1^1+{\psi }_1^k-{\varphi }_2^2)& i\in S_2^k\\ f({\varphi }_i^i-{\psi }_{i-1}^k,{\varphi }_i^i+{\psi }_i^k-2\pi )& i\in S_3^k\\ f({\varphi }_i^i-{\psi }_{i-1}^k+2\pi -{\varphi }_{i-1}^{i-1},{\varphi }_i^i+{\psi }_i^k-{\varphi }_{i+1}^{i+1})& i\in S_4^k\\ f({\varphi }_i^i-{\psi }_{i-1}^k,{\varphi }_i^i+{\psi }_i^k-{\varphi }_{i+1}^{i+1})& i\in S_5^k\\ f({\varphi }_i^i-{\psi }_{i-1}^k+2\pi -{\varphi }_{i-1}^{i-1},{\varphi }_i^i+{\psi }_i^k-2\pi )& i\in S_6^k\\ f_{\text{B}}({\varphi }_{N+1}^{N+1}-{\psi }_N^k)& i\in S_7^k\\ f_{\text{B}}({\varphi }_{N+1}^{N+1}-{\psi }_N^k+2\pi -{\varphi }_N^N)& i\in S_8^k\end{array}& (14)\end{array}$

where $S_{1-8}^k$ 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, $S_{1,2,7,8}^k$ are for i = 1 and N + 1, which are given by

$\begin{array}{l}{S}_1^k=\left\{i\mid i=1\right\}\cap \left\{i\mid {\varphi }_i^k\in ({\varphi }_i^{i+1},{\varphi }_i^i]\right\}\\ S_2^k=\left\{i\mid i=1\right\}\cap \left\{i\mid i\notin S_1^k\right\}\\ S_7^k=\left\{i\mid i=N+1\right\}\cap \left\{i\mid {\varphi }_i^k\in ({\varphi }_i^{i-1},{\varphi }_i^i]\right\}\\ S_8^k=\left\{i\mid i=N+1\right\}\cap \left\{i\mid i\notin S_7^k\right\},\end{array}$

$S_{3,4,5,6}^k$ are for i = 2 to N, which are given by

$\begin{array}{l}{S}_3^k=\left\{i\mid i=k\right\}\cup \left\{i\in \left[2,k-1\right]\mid {\varphi }_i^k\in ({\varphi }_i^{i+1},{\varphi }_i^i]\right\}\\ \cup \left\{i\in \left[k+1,N\right]\mid {\varphi }_i^k\in ({\varphi }_i^{i-1},{\varphi }_i^i]\right\}\\ S_4^k=\left\{i\in \left[2,k-1\right]\mid {\varphi }_{i-1}^k\in ({\varphi }_{i-1}^i,{\varphi }_{i-1}^{i-1}]\right\}\\ \cup \left\{i\in \left[k+1,N\right]\mid {\varphi }_{i+1}^k\in ({\varphi }_{i+1}^i,{\varphi }_{i+1}^{i+1}]\right\}\\ S_5^k=\left\{i\in \left[2,k-1\right]\mid i\notin S_3^k\text{ and }i\notin S_4^k\right\}\\ S_6^k=\left\{i\in \left[k+1,N\right]\mid i\notin S_3^k\text{ and }i\notin S_4^k\right\}.\end{array}$

From (11), periodic solutions must satisfy

$\begin{array}{ll}{\widehat{\varphi }}_1^1={\widehat{\varphi }}_2^2=\cdots ={\widehat{\varphi }}_{N+1}^{N+1}\equiv {\widehat{\varphi }}^{\text{td}}.& (15)\end{array}$

This means that the touchdown phase of all oscillators equals ${\widehat{\varphi }}^{\text{td}}$. From (14) and (15), we obtain the periodic solution by ${\widehat{\varphi }}^{\text{td}}$ and ${\widehat{\psi }}_i^k(i=1,\dots ,N)$ as follows:

$\begin{array}{ll}{\widehat{\varphi }}^{\text{td}}=2\pi -m_{\text{B}}b& (16)\end{array}$

$\begin{array}{ll}{\widehat{\psi }}_i^k=\{\begin{array}{cc}{\widehat{\psi }}_{k-1}^k& i=k-1\\ 2(1-\beta )\pi -(2m_{\text{B}}-m)b& i\in T_1^k\\ 2(1-\beta )\pi -(3m_{\text{B}}-m)b& i\in \left[1,k-2\right]\text{ and }i\notin T_1^k\\ {\widehat{\psi }}_k^k& i=k\\ 2\beta \pi +(2m_{\text{B}}-m)b& i+1\in T_2^k\\ 2\beta \pi +(3m_{\text{B}}-m)b& \text{otherwise},\end{array}& (17)\end{array}$

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 $T_{1,2}^k$ are given by

$\begin{array}{l}{T}_1^k=\{i\in \left[1,k-2\right]|2\beta \pi +(2m_{\text{B}}-m)b<{\widehat{\varphi }}_i^k\le {\widehat{\varphi }}^{\text{td}}\}\\ T_2^k=\{i\in \left[k+2,N+1\right]|2\beta \pi +(2m_{\text{B}}-m)b<{\widehat{\varphi }}_i^k\le {\widehat{\varphi }}^{\text{td}}\}.& (18)\end{array}$

$T_1^k$ 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., $S_3^k$ in the direct-wave region). $T_2^k$ 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., $S_3^k$ in the retrograde-wave region). ${\widehat{\psi }}_{k-1}^k$ is not determined uniquely but satisfies

$\begin{array}{ll}2(1-\beta )\pi -(3m_{\text{B}}-m)b<{\widehat{\psi }}_{k-1}^k\le 2\beta \pi .& (19)\end{array}$

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). ${\widehat{\psi }}_k^k$ is represented using ${\widehat{\psi }}_{k-1}^k$ as

$\begin{array}{ll}{\widehat{\psi }}_k^k=\{\begin{array}{cc}{\widehat{\psi }}_{k-1}^k+(4\beta -2)\pi +(4m_{\text{B}}-m)b& {\psi }_{\text{A}}<{\widehat{\psi }}_{k-1}^k<{\psi }_{\text{B}}\\ 2\beta \pi +(3m_{\text{B}}-m)b& {\psi }_{\text{B}}\le {\widehat{\psi }}_{k-1}^k\le 2\beta \pi \end{array}& (20)\end{array}$

where ψ_A = 2(1 − β)π − (3m_B − m)b and ψ_B = 2(1 − β)π − m_Bb. Therefore, the relationship of ${\widehat{\psi }}_{k-1}^k$ and ${\widehat{\psi }}_k^k$ 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 $S_{1-8}^k$.

FIGURE 10

Figure 10. Relation between ${\widehat{\psi }}_{k-1}^k$ and ${\widehat{\psi }}_k^k$.

3.4. Stability of Source-Wave Gait

To discuss stability, we added perturbation $\Delta {\psi }_i^k(i=1,\dots ,N)$ to the obtained periodic solutions on Σ^k and evaluated the linear stability by analytically calculating the linear map of perturbation:

$\begin{array}{ll}\Delta {\Psi }^k\mapsto P_k\Delta {\Psi }^k,& (21)\end{array}$

where $\Delta {\Psi }^k={\left[\Delta {\psi }_1^k\Delta {\psi }_2^k\dots \Delta {\psi }_N^k\right]}^T$, and $P_k\in {ℝ}^{N\times 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 ${\widehat{\varphi }}^{\mathrm{\text{td}}}+\Delta {\varphi }_i^i$ and substitute ${\psi }_i^k={\widehat{\psi }}_i^k+\Delta {\psi }_i^k$ and ${\varphi }_i^i={\widehat{\varphi }}^{\text{td}}+\Delta {\varphi }_i^i$ into (14), $\Delta {\varphi }_i^i$ (i = 1, …, N + 1) can be represented with $\Delta {\psi }_i^k$. The matrix P_k is derived by substituting ${\varphi }_i^i={\widehat{\varphi }}^{\text{td}}+\Delta {\varphi }_i^i$ into (11). The eigenvalues λ_Sk of P_k for the source-wave gait are derived depending on solution sets A and B by

$\begin{array}{ll}{\lambda }_{\text{S}k}=\{\begin{array}{cc}\left\{1,1/2,2/3,\dots ,2/3\right\}& {\psi }_{\text{A}}<{\widehat{\psi }}_{k-1}^k<{\psi }_{\text{B}}\\ \left\{1,2/3,2/3,\dots ,2/3\right\}& {\psi }_{\text{B}}\le {\widehat{\psi }}_{k-1}^k\le 2\beta \pi \end{array}& (22)\end{array}$

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 ${\varphi }_{k-1}^k$ 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 $({\tilde{\varphi }}_1^k,{\tilde{\varphi }}_2^k,\dots ,{\tilde{\varphi }}_{N+1}^k)$ just before the touchdown of Leg k by ${\Phi }_i(t;{\tilde{\varphi }}_1^k,{\tilde{\varphi }}_2^k,\dots ,{\tilde{\varphi }}_{N+1}^k)$. From the front–rear symmetry of our model, the following equation is satisfied for the initial phases such that ${\tilde{\varphi }}_i^k={\tilde{\varphi }}_i^{N+2-k}$ for i = 1, …, N + 1:

$\begin{array}{l}{\Phi }_i(t;{\tilde{\varphi }}_1^k,{\tilde{\varphi }}_2^k,\dots ,{\tilde{\varphi }}_{N+1}^k)\\ ={\Phi }_{N+2-i}(t;{\tilde{\varphi }}_{N+1}^{N+2-k},{\tilde{\varphi }}_N^{N+2-k},\dots ,{\tilde{\varphi }}_1^{N+2-k})& (23)\end{array}$

The periodic solution of the direct-wave gait ${\widehat{\psi }}_i^N(i=1,\dots ,N)$ on Σ^N, where ${\widehat{\varphi }}^{\text{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 ${\Phi }_i(t;{\tilde{\varphi }}_{N+1}^N,\dots ,{\tilde{\varphi }}_1^N)={\Phi }_{N+2-i}(t;{\tilde{\varphi }}_1^2,\dots ,{\tilde{\varphi }}_{N+1}^2)$. When the periodic solution ${\widehat{\varphi }}_j^N$ (j = 1, …, N + 1) is used for the initial value ${\tilde{\varphi }}_j^N$, we obtain ${\widehat{\varphi }}_i^N={\widehat{\varphi }}_{N+2-i}^2$ (i = 1, …, N). This yields ${\widehat{\psi }}_i^N=2\pi -{\widehat{\psi }}_{N+1-i}^2$ because ${\widehat{\psi }}_i^N={\widehat{\varphi }}_{i+1}^N-{\widehat{\varphi }}_i^N$. By substituting solution set B of the retrograde-wave gait ${\widehat{\psi }}_{N+1-i}^2$ (k = 2 in (17)), we obtain the solution of the direct-wave gait as

$\begin{array}{ll}{\widehat{\psi }}_i^N=\{\begin{array}{cc}{\widehat{\psi }}_N^N& i=N\\ 2(1-\beta )\pi -(2m_{\text{B}}-m)b& i\in T_1^N\\ 2(1-\beta )\pi -(3m_{\text{B}}-m)b& \text{otherwise}.\end{array}& (24)\end{array}$

where ${\widehat{\psi }}_N^N$ is an arbitrary constant fulfilling

$\begin{array}{ll}2(1-\beta )\pi \le {\widehat{\psi }}_N^N\le 2\beta \pi +m_{\text{B}}b.& (25)\end{array}$

The linear stability of the solution can also be calculated using the symmetry condition (23). The symmetry condition gives the relation $\Delta {\psi }_i^N=-\Delta {\psi }_{N+1-i}^2$, 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 $P_{\text{D}}=Q{P}_2{Q}^{-1}$. Thus, the eigenvalues λ_D of P_D are obtained by those of P₂ as

$\begin{array}{ll}{\lambda }_{\text{D}}=\left\{1,2/3,2/3,\dots ,2/3\right\}.& (26)\end{array}$

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 $T_1^k$, 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 ${\widehat{\psi }}_s^k$ 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 ${\widehat{\psi }}_{i-1}^i$-${\widehat{\psi }}_i^i$-${\widehat{\psi }}_{i+1}^i$ space; (B) direct- and source-wave gaits (k = N−1, N) projected onto the ${\widehat{\psi }}_{N-2}^{N-1}$-${\widehat{\psi }}_{N-1}^{N-1}$-${\widehat{\psi }}_N^{N-1}$ space; and (C) retrograde- (k = 2) and source-wave gaits (k = 3, 4) projected onto the ${\widehat{\psi }}_1^2$-${\widehat{\psi }}_2^2$-${\widehat{\psi }}_3^2$ 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

$\begin{array}{ll}\xi =\frac{\frac{1}{1-\beta }-\frac{m_{\text{B}}g}{2\kappa a}}{1-\frac{(3m_{\text{B}}-m)g}{2\kappa a}}\eta ,& (27)\end{array}$

$\begin{array}{ll}v=\frac{\xi \omega }{2\pi }.& (28)\end{array}$

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

When k = 2 and ${\widehat{\psi }}_1\approx 2\beta \pi$, 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 ${\widehat{\psi }}_N\approx 2(1-\beta )\pi$ 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 ($\sqrt{\kappa /m_i}\gg \omega$). 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

$\begin{array}{ll}{E}^j={\displaystyle \sum _{i=1}^N}{\psi }_i^j.& (29)\end{array}$

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

$\begin{array}{ll}{E}^j\mapsto E^j+{\varphi }_1^1-{\varphi }_{N+1}^{N+1}.& (30)\end{array}$

The relation ${\varphi }_1^1={\varphi }_{N+1}^{N+1}={\widehat{\varphi }}^{\text{td}}$ holds if the neighboring legs (Legs 2 and N) are in the stance phases (${\varphi }_2^1\in \left[0,2\beta \pi \right]$ and ${\varphi }_N^{N+1}\in \left[0,2\beta \pi \right]$) 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 ${\varphi }_1^1$ and ${\varphi }_{N+1}^{N+1}$ using (3) and (12) when ${\varphi }_2^1\in \left[0,2\beta \pi \right]$ and ${\varphi }_N^{N+1}\in \left[0,2\beta \pi \right]$.

$\begin{array}{ll}{\varphi }_1^1=2\pi -\frac{(1-\beta )m_{1}g\pi }{\kappa a},& (31)\end{array}$

$\begin{array}{ll}{\varphi }_{N+1}^{N+1}=2\pi -\frac{(1-\beta )m_{N+1}g\pi }{\kappa a}.& (32)\end{array}$

These equations show that ${\varphi }_1^1={\varphi }_{N+1}^{N+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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar