††thanks: Corresponding author. E-mail<EMAIL_ADDRESS>
# A Central Pattern Generator Network for Simple Control of Gait Transitions
in Hexapod Robots based on Phase Reduction

Norihisa Namura  Department of Systems and Control Engineering, Tokyo
Institute of Technology, Tokyo 152-8552, Japan  Hiroya Nakao  Department of
Systems and Control Engineering, Tokyo Institute of Technology, Tokyo
152-8552, Japan  Research Center for Autonomous Systems Materialogy, Institute
of Innovative Research, Tokyo Institute of Technology, Yokohama, 226-8501,
Japan

###### Abstract

We present a model of the central pattern generator (CPG) network that can
control gait transitions in hexapod robots in a simple manner based on phase
reduction. The CPG network consists of six weakly coupled limit-cycle
oscillators, whose synchronization dynamics can be described by six phase
equations through phase reduction. Focusing on the transitions between the
hexapod gaits with specific symmetries, the six phase equations of the CPG
network can further be reduced to two independent equations for the phase
differences. By choosing appropriate coupling functions for the network, we
can achieve desired synchronization dynamics regardless of the detailed
properties of the limit-cycle oscillators used for the CPG. The effectiveness
of our CPG network is demonstrated by numerical simulations of gait
transitions between the wave, tetrapod, and tripod gaits, using the FitzHugh-
Nagumo oscillator as the CPG unit.

##  I Introduction

Legged locomotion is a key ability for many terrestrial organisms. For
instance, humans use bipedal locomotion, cats and horses use quadruped
locomotion, and insects use hexapodal locomotion. To adapt to various terrains
or to change the speed of locomotion, gait transitions often occur in
vertebrates (e.g., horses Heglund1988speed ) and invertebrates (e.g.,
drosophilas Strauss1990coordination  and stick insects Graham1981coordinated
).

Some neurophysiological experiments, such as the one for decerebrated cats
Shik1966control , suggest that gait patterns are generated in the neuronal
network called the central pattern generator (CPG) in the spinal cord
Grillner1975locomotion . In invertebrates, CPGs also exist as distributed
neuronal circuits that control gait patterns Marder2005invertebrate ;
Grillner2006biological . The CPG can produce self-sustained rhythms for
generating basic gait patterns without external stimulation or sensory
feedback, and it can further be modulated by the signals from the brain or
sensory organs to switch between various rhythmic patterns, i.e., gait
patterns Rybak2006modelling .

Compared to wheeled robots, legged robots have the advantage of adapting to
uneven terrains and physical environments more efficiently, and many studies
have been conducted on biped robots Aoi2005locomotion , quadruped robots
Owaki2017quadruped , hexapod robots Ambe2013stability ; Minati2018versatile ,
and salamander robots Ijspeert2007swimming . In particular, hexapod robots are
more stable than biped or quadruped robots and have been applied, e.g., to
humanitarian demining Hexapod_demine  and exploration of the surface of Mars
Hexapod_mars . Since gait transitions are also important for hexapod robots,
methods for realizing typical insect gaits and transitions between them using
hexapod robots have been studied Inagaki2006wave ; Ricardo2010hexapod ;
Chen2012smooth ; Manoonpong2013neural ; Yu2016gait ; Bal2021neural ;
Zhu2022generic .

For controlling the locomotion of legged robots, CPG-based control strategy is
an efficient bio-inspired method because it provides self-oscillations and
requires less computational costs Bal2021neural . CPG networks are widely
modeled as coupled limit-cycle oscillators. In rhythmic systems of mutually
coupled limit-cycle oscillators like the CPG networks, synchronization can
occur, leading to the alignment of their rhythms with each other.
Synchronization of weakly coupled limit-cycle oscillators can be analyzed by
the phase reduction theory Kuramoto1984chemical ; Hoppensteadt1997weakly ;
Winfree2001geometry ; Ermentrout2010mathematical ; Nakao2016phase , which can
approximately reduce the nonlinear multidimensional dynamics of the oscillator
to a single-variable phase equation characterized by the natural frequency and
phase response property of the oscillator. Phase reduction can also be useful
for engineering limit-cycle oscillators, e.g., in optimizing entrainment
Harada2010optimal ; Dasanayake2011optimal ; Zlotnik2013optimal ;
Monga2019phase ; Kato2021optimization ; Takata2021fast  and mutual
synchronization Shirasaka2017optimizing ; Watanabe2019optimization ;
Namura2024optimal . Moreover, it has been used for analyzing gait-generation
models Ghigliazza2004minmal  and gait transitions Aminzare2018gait ;
Yeldesbay2018role  in hexapods.

In this study, we present a CPG network for controlling gait transitions in
hexapod robots, which allows various types of smooth limit-cycle oscillators
to be used as the CPG unit. We use phase reduction to derive six phase
equations from the mathematical model of the CPG network consisting of six
limit-cycle oscillators with weak coupling. Using the symmetries of the
typical hexapod gaits, i.e., the wave, tetrapod, and tripod gaits, we further
reduce the six phase equations to two independent phase-difference equations.
By choosing appropriate coupling functions for the network, we can achieve all
typical gaits irrespective of the detailed properties of the limit-cycle
oscillators. Compared to existing studies, our CPG network provides a simpler
and more general framework in which the CPG model and coupling functions can
be optimized. We demonstrate the effectiveness of our CPG network by numerical
simulations of gait transitions between the wave, tetrapod, and tripod gaits,
adopting the FitzHugh-Nagumo oscillator as the CPG unit.

This paper is organized as follows. We first present a CPG network and explain
the symmetric typical hexapod gaits in Sec. II. We next describe the gait
control methods in Sec. III. We demonstrate the results of the gait
transitions by the present CPG network in Sec. IV, and we give conclusions in
Sec. V. Appendix A gives further examples of the symmetric gait patterns.

##  II CPG Network

In this section, we present a ladder-type bidirectional CPG network consisting
of six limit-cycle oscillators as schematically shown in Fig. 1.

Figure 1:  Coupling scheme of the present CPG network. The six CPGs are
denoted by 1 to 6, which correspond to LF, LM, LH, RF, RM, and RH,
respectively. The weights $c_{1,2}$ determine the coupling strength and
$b_{1,2}$ determine whether the coupling is excitatory or inhibitory. The
functions $\Gamma_{1,2}$ determine the phase coupling between the CPGs. The
red solid arrows show the coupling by $\Gamma_{1}$ and the blue dashed arrows
show the coupling by $\Gamma_{2}$, respectively.

###  II.1 Phase dynamics of the CPG network

We assume that the dynamics of the six limit-cycle oscillators in the CPG
network (Fig. 1) are described by

$\displaystyle\begin{aligned}
\dot{\bm{X}}_{1}={}&s\bm{F}(\bm{X}_{1})+\varepsilon\left(b_{1}c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{1},\bm{X}_{4})+b_{2}c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{1},\bm{X}_{2})\right),\\\
\dot{\bm{X}}_{2}={}&s\bm{F}(\bm{X}_{2})+\varepsilon
b_{1}c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{2},\bm{X}_{5})\\\
&+\varepsilon\left(c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{2},\bm{X}_{1})+b_{2}c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{2},\bm{X}_{3})\right),\\\
\dot{\bm{X}}_{3}={}&s\bm{F}(\bm{X}_{3})+\varepsilon\left(b_{1}c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{3},\bm{X}_{6})+c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{3},\bm{X}_{2})\right),\\\
\dot{\bm{X}}_{4}={}&s\bm{F}(\bm{X}_{4})+\varepsilon\left(c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{4},\bm{X}_{1})+b_{2}c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{4},\bm{X}_{5})\right),\\\
\dot{\bm{X}}_{5}={}&s\bm{F}(\bm{X}_{5})+\varepsilon
c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{5},\bm{X}_{2})\\\
&+\varepsilon\left(c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{5},\bm{X}_{4})+b_{2}c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{5},\bm{X}_{6})\right),\\\
\dot{\bm{X}}_{6}={}&s\bm{F}(\bm{X}_{6})+\varepsilon\left(c_{1}\tilde{\bm{H}}_{1}(\bm{X}_{6},\bm{X}_{3})+c_{2}\tilde{\bm{H}}_{2}(\bm{X}_{6},\bm{X}_{5})\right),\end{aligned}$
(1)

where $\bm{X}_{k}(t)\in\mathbb{R}^{N}\;(k=1,2,3,4,5,6)$ denotes each
oscillator state at time $t$, $\dot{\bm{X}}$ represents the time derivative of
$\bm{X}$, $\bm{F}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is a smooth vector field
common to all oscillators, $s>0$ is a common timescale parameter that
determines the oscillation frequency of all the oscillators, $c_{1,2}>0$ are
the coupling strength of $O(1)$, $b_{1,2}\in\\{-1,+1\\}$ determine whether the
coupling is excitatory or inhibitory, $0<\varepsilon\ll 1$ is a small
parameter representing weakness of the coupling, and
$\tilde{\bm{H}}_{1,2}:\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$
represent smooth coupling functions between the oscillators. We assume that
the original system $\dot{\bm{X}}=\bm{F}(\bm{X})$ without scaling ($s=1$) has
an exponentially stable limit cycle $\chi$ of natural period $T$ and frequency
$\omega=2\pi/T$, where the state on $\chi$ at time $t$ is represented as
$\tilde{\bm{X}}_{0}(t)$ and satisfies
$\tilde{\bm{X}}_{0}(t)=\tilde{\bm{X}}_{0}(t+T)$, given the origin state
$\tilde{\bm{X}}_{0}(0)$ on the limit cycle. Then, the period and frequency of
the limit cycle $\chi$ of the time-rescaled vector field $s\bm{F}$ are $T/s$
and $s\omega$, respectively.

The mathematical model of weakly coupled identical limit-cycle oscillators can
be analyzed by phase reduction theory Kuramoto1984chemical ;
Hoppensteadt1997weakly ; Winfree2001geometry ; Ermentrout2010mathematical ;
Nakao2016phase . We first introduce an asymptotic phase function
$\Theta:\mathcal{B}\to[0,2\pi)$ into the basin
$\mathcal{B}\subset\mathbb{R}^{N}$ of $\chi$ satisfying
$\nabla\Theta(\bm{X})\cdot\bm{F}(\bm{X})=\omega$ for all
$\bm{X}\in\mathcal{B}$, where $\bm{a}\cdot\bm{b}=\sum_{n=1}^{N}a_{n}b_{n}$
represents the scalar product of two vectors $\bm{a}$,
$\bm{b}\in\mathbb{R}^{N}$. Using the phase function, we can define the phase
value $\theta$ of the oscillator state $\bm{X}\in\mathcal{B}$ by
$\theta=\Theta(\bm{X})$, which increases constantly at the natural frequency
$\omega$ over time,

$\displaystyle\dot{\theta}(t)=\nabla\Theta(\bm{X}(t))\cdot\bm{F}(\bm{X}(t))=\omega,$
(2)

where we regard that the phase values $0$ and $2\pi$ are identical. The state
on $\chi$ can be represented by
$\bm{X}_{0}(\theta)=\tilde{\bm{X}}_{0}(t=\theta/\omega)$ as a $2\pi$-periodic
function of $\theta$, which satisfies
$\bm{X}_{0}(\theta)=\bm{X}_{0}(\theta+2\pi)$.

Since all coupling terms are of $O(\varepsilon)$ and sufficiently weak, we can
assume that the state of each limit-cycle oscillator is close to $\chi$, i.e.,

$\displaystyle\bm{X}_{k}(t)=\bm{X}_{0}(\theta_{k}(t))+O(\varepsilon)\quad(k=1,2,3,4,5,6),$
(3)

where $\theta_{k}=\Theta(\bm{X}_{k})$ for $k=1,2,3,4,5,6$. Substituting Eq.
(3) into Eq. (1) and neglecting errors of $O(\varepsilon^{2})$, we can derive
the following phase equations:

$\displaystyle\begin{aligned}
\dot{\theta}_{1}={}&s\omega+\varepsilon\bm{Z}(\theta_{1})\cdot\left(b_{1}c_{1}\bm{H}_{1}(\theta_{1},\theta_{4})+b_{2}c_{2}\bm{H}_{2}(\theta_{1},\theta_{2})\right),\\\
\dot{\theta}_{2}={}&s\omega+\varepsilon\bm{Z}(\theta_{2})\cdot
b_{1}c_{1}\bm{H}_{1}(\theta_{2},\theta_{5}),\\\
&+\varepsilon\bm{Z}(\theta_{2})\cdot\left(c_{2}\bm{H}_{2}(\theta_{2},\theta_{1})+b_{2}c_{2}\bm{H}_{2}(\theta_{2},\theta_{3})\right),\\\
\dot{\theta}_{3}={}&s\omega+\varepsilon\bm{Z}(\theta_{3})\cdot\left(b_{1}c_{1}\bm{H}_{1}(\theta_{3},\theta_{6})+c_{2}\bm{H}_{2}(\theta_{3},\theta_{2})\right),\\\
\dot{\theta}_{4}={}&s\omega+\varepsilon\bm{Z}(\theta_{4})\cdot\left(c_{1}\bm{H}_{1}(\theta_{4},\theta_{1})+b_{2}c_{2}\bm{H}_{2}(\theta_{4},\theta_{5})\right),\\\
\dot{\theta}_{5}={}&s\omega+\varepsilon\bm{Z}(\theta_{5})\cdot
c_{1}\bm{H}_{1}(\theta_{5},\theta_{2}),\\\
&+\varepsilon\bm{Z}(\theta_{5})\cdot\left(c_{2}\bm{H}_{2}(\theta_{5},\theta_{4})+b_{2}c_{2}\bm{H}_{2}(\theta_{5},\theta_{6})\right),\\\
\dot{\theta}_{6}={}&s\omega+\varepsilon\bm{Z}(\theta_{6})\cdot\left(c_{1}\bm{H}_{1}(\theta_{6},\theta_{3})+c_{2}\bm{H}_{2}(\theta_{6},\theta_{5})\right).\end{aligned}$
(4)

Note that the frequency of each oscillator is scaled to $s\omega$. Here,
$\bm{Z}:[0,2\pi)\to\mathbb{R}^{N}$ is the phase sensitivity function (PSF) of
the oscillator defined by
$\bm{Z}(\theta)=\left.\nabla\Theta(\bm{X})\right|_{\bm{X}=\bm{X}_{0}(\theta)}$,
which characterizes the linear response of the oscillator phase $\theta$ to
small perturbations. The PSF can be calculated as the $2\pi$-periodic solution
to the following adjoint equation Ermentrout2010mathematical ; Brown2004phase
:

$\displaystyle\omega\frac{d}{d\theta}\bm{Z}(\theta)=-\bm{J}(\bm{X}_{0}(\theta))^{\top}\bm{Z}(\theta),$
(5)

where $\bm{J}:\mathbb{R}^{N}\to\mathbb{R}^{N\times N}$ is the Jacobian matrix
of the vector field $\bm{F}$, i.e., $\bm{J}(\bm{X})=\nabla\bm{F}(\bm{X})$. The
PSF satisfies the following normalization condition:

$\displaystyle\bm{Z}(\theta)\cdot\frac{d\bm{X}_{0}(\theta)}{d\theta}=1.$ (6)

The functions $\bm{H}_{1,2}:[0,2\pi)\times[0,2\pi)\to\mathbb{R}^{N}$ are
mutual coupling functions (MCFs) defined by

$\displaystyle\bm{H}_{1,2}(\theta_{i},\theta_{j})$
$\displaystyle=\tilde{\bm{H}}_{1,2}(\bm{X}_{0}(\theta_{i}),\bm{X}_{0}(\theta_{j})),$
(7)

which depend on two phases $\theta_{i}$ and $\theta_{j}$ for
$i,j\in\\{1,2,3,4,5,6\\}$. From Eq. (3),

$\displaystyle\tilde{\bm{H}}_{1,2}(\bm{X}_{i},\bm{X}_{j})$
$\displaystyle=\bm{H}_{1,2}(\theta_{i},\theta_{j})+O(\varepsilon)$ (8)

holds for $i,j\in\\{1,2,3,4,5,6\\}$.

Since the frequency of the oscillator is $s\omega$ and therefore the relative
phase $\theta_{k}-s\omega t$ is a slow variable, we can conduct the averaging
approximation Kuramoto1984chemical ; Hoppensteadt1997weakly ; Nakao2016phase
over one period $T/s$ and obtain the following averaged phase equations by
neglecting errors of $O(\varepsilon^{2})$:

$\displaystyle\begin{aligned}
\dot{\theta}_{1}={}&s\omega+\varepsilon(b_{1}c_{1}\Gamma_{1}(\theta_{1}-\theta_{4})+b_{2}c_{2}\Gamma_{2}(\theta_{1}-\theta_{2})),\\\
\dot{\theta}_{2}={}&s\omega+\varepsilon
b_{1}c_{1}\Gamma_{1}(\theta_{2}-\theta_{5})\\\
&+\varepsilon(c_{2}\Gamma_{2}(\theta_{2}-\theta_{1})+b_{2}c_{2}\Gamma_{2}(\theta_{2}-\theta_{3})),\\\
\dot{\theta}_{3}={}&s\omega+\varepsilon(b_{1}c_{1}\Gamma_{1}(\theta_{3}-\theta_{6})+c_{2}\Gamma_{2}(\theta_{3}-\theta_{2})),\\\
\dot{\theta}_{4}={}&s\omega+\varepsilon(c_{1}\Gamma_{1}(\theta_{4}-\theta_{1})+b_{2}c_{2}\Gamma_{2}(\theta_{4}-\theta_{5})),\\\
\dot{\theta}_{5}={}&s\omega+\varepsilon
c_{1}\Gamma_{1}(\theta_{5}-\theta_{2})\\\
&+\varepsilon(c_{2}\Gamma_{2}(\theta_{5}-\theta_{4})+b_{2}c_{2}\Gamma_{2}(\theta_{5}-\theta_{6})),\\\
\dot{\theta}_{6}={}&s\omega+\varepsilon(c_{1}\Gamma_{1}(\theta_{6}-\theta_{3})+c_{2}\Gamma_{2}(\theta_{6}-\theta_{5})),\\\
\end{aligned}$ (9)

where the functions $\Gamma_{1,2}$ are called phase coupling functions (PCFs),
which are $2\pi$-periodic functions represented by

$\displaystyle\Gamma_{1,2}(\varphi)=\left\langle\bm{Z}(\psi)\cdot\bm{H}_{1,2}(\psi,\psi-\varphi)\right\rangle_{\psi}.$
(10)

Here, we denote the averaging of a smooth $2\pi$-periodic function $g$ over
one period of oscillation as

$\displaystyle\left\langle
g(\psi)\right\rangle_{\psi}=\frac{1}{2\pi}\int_{0}^{2\pi}g(\psi)d\psi.$ (11)

Now, we assume that the MCFs are given in the following specific form:

$\displaystyle\bm{H}_{1,2}(\theta_{i},\theta_{j})=\frac{\bm{Z}(\theta_{i})P_{1,2}(\theta_{i}-\theta_{j})}{\left\langle\left\|\bm{Z}(\psi)\right\|^{2}\right\rangle_{\psi}},$
(12)

where $P_{1,2}:\mathbb{R}\to\mathbb{R}$ are arbitrary $2\pi$-periodic smooth
functions. As discussed in Namura2024optimal , these MCFs give the most
energy-efficient forms of the mutual coupling that maximize the absolute value
of the corresponding PCFs at each value of the phase difference under a given
power, leading to fast mutual synchronization of the oscillators. Introducing
$\varphi=\theta_{i}-\theta_{j}$, the PCFs are obtained from the above MCFs as

$\displaystyle\Gamma_{1,2}(\varphi)=P_{1,2}(\varphi).$ (13)

By choosing the functional forms of $P_{1,2}(\varphi)$, we can realize any
forms of $\Gamma_{1,2}$. Thus, we need to consider only the PCFs
$\Gamma_{1,2}$ in designing the CPG network, and then MCFs $\bm{H}_{1,2}$
follow from Eq. (12). We note that, in general, the functional form of the PCF
depends on both the PSF and MCF. For the above particular MCF, the PCF become
independent of the PSF and we can arbitrarily design it to realize desired
synchronization dynamics.

###  II.2 Symmetric gait patterns

Next, we describe the hexapod gaits considered in this study. The six legs of
the hexapod robots are called the left front (LF), left middle (LM), left hind
(LH), right front (RF), right middle (RM), and right hind (RH), as shown in
Fig. 1. We consider
$(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6})$ to be
the phases of (LF, LM, LH, RF, RM, RH), respectively.

We assume that the hexapod gaits have the following symmetries:

$\displaystyle\theta_{6}-\theta_{3}=\theta_{5}-\theta_{2}=\theta_{4}-\theta_{1}$
$\displaystyle\mod 2\pi,$ (14)
$\displaystyle\theta_{3}-\theta_{2}=\theta_{2}-\theta_{1}=\theta_{6}-\theta_{5}=\theta_{5}-\theta_{4}$
$\displaystyle\mod 2\pi,\begin{aligned} \end{aligned}$ (15)

where Eq. (14) indicates that the phase differences between the contralateral
legs in the front, middle, and hind positions are equal, and Eq. (15)
indicates that the phase differences between the adjacent ipsilateral legs are
equal, i.e., the phase shift of the middle legs from the front legs is the
same as that of the hind legs from the middle legs on both sides.

As an example, we consider three typical hexapod gaits, i.e., the slowest gait
called the wave gait, the middle-speed gait called the tetrapod gait, and the
fastest gait called the tripod gait Wilson1966insect ; Collins1993hexapodal ;
Golubitsky1998modular . Figure 2 shows the diagrams of the three gaits, where
the states (stance or swing) of the six legs are plotted in gray and black
colors, respectively. In the wave gait (also called a metachronal gait), the
hexapod robot moves the six legs, i.e., (RH), (RM), (RF), (LH), (LM), and
(LF), in order at equal phase intervals of $\pi/3$. In the tetrapod gait, the
legs form three synchronized pairs, i.e., (LH, RM), (LM, RF), and (LF, RH),
which are moved in order at equal phase intervals of $2\pi/3$. Finally, in the
tripod gait, the legs form two synchronized triplets, (RF, LM, RH) and (LF,
RM, LH), which are alternatively moved (i.e., with the phase interval of
$\pi$).

The phase shifts of the legs in the three types of the gaits correspond to

$\displaystyle\begin{aligned}
&(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6})\\\
&=\left\\{\begin{array}[]{ll}\left(0,\frac{1}{3}\pi,\frac{2}{3}\pi,\pi,\frac{4}{3}\pi,\frac{5}{3}\pi\right)\quad(\textrm{wave}),\\\
\\\
\left(0,\frac{2}{3}\pi,\frac{4}{3}\pi,\frac{2}{3}\pi,\frac{4}{3}\pi,0\right)\quad(\textrm{tetrapod}),\\\
\\\
\left(0,\pi,0,\pi,0,\pi\right)\quad(\textrm{tripod}),\end{array}\right.\end{aligned}$
(16)

when the LF is chosen as the reference, $\theta_{1}=0$. We can find that all
the three gaits satisfy the symmetries given in Eqs. (14) and (15). We note
that other hexapod gaits are also possible as described in Appendix A.

Figure 2: Diagrams of the wave, tetrapod, and tripod gaits (see Sec. III A for
the definition of the swing and stance phases).

###  II.3 Equations for the phase differences

We now show that the six phase equations can be reduced to two independent
phase-difference equations by using the symmetries (14) and (15) of the gaits.
First, we consider the symmetry of the phase differences between the
contralateral legs given in Eq. (14), which yields

$\displaystyle\theta_{k+3}-\theta_{m+3}=\theta_{k}-\theta_{m}\mod
2\pi\quad(k,m=1,2,3).$ (17)

Plugging Eqs. (14) and (17) into Eq. (9), we find that the time derivatives of
$\theta_{k+3}-\theta_{k}$ take the same value for $k=1,2,3$, i.e.,

$\displaystyle\dot{\theta}_{k+3}-\dot{\theta}_{k}=\varepsilon
c_{1}\left(\Gamma_{1}(\theta_{k+3}-\theta_{k})-b_{1}\Gamma_{1}(\theta_{k}-\theta_{k+3})\right).$
(18)

Therefore, if the initial phases $\\{\theta_{k}(0)\\}_{k=1}^{6}$ satisfy the
condition (14), the phases satisfy

$\displaystyle\theta_{k+3}(t)-\theta_{k}(t)\coloneqq\alpha(t)\mod
2\pi\quad(k=1,2,3)$ (19)

for all $t$, i.e., the phase differences $\theta_{4}-\theta_{1}$,
$\theta_{5}-\theta_{2}$, and $\theta_{6}-\theta_{3}$ are always equal to
$\alpha$. Therefore, the time evolution of the phase difference $\alpha$ obeys

$\displaystyle\dot{\alpha}=\varepsilon
c_{1}\left(\Gamma_{1}(\alpha)-b_{1}\Gamma_{1}(-\alpha)\right).$ (20)

From Eqs. (17) and (19), by using the phase difference $\alpha$, the dynamics
of the six phases can be rewritten as

$\displaystyle\begin{aligned}
\dot{\theta}_{1}&=s\omega+\varepsilon(b_{1}c_{1}\Gamma_{1}(-\alpha)+b_{2}c_{2}\Gamma_{2}(\theta_{1}-\theta_{2})),\\\
\dot{\theta}_{2}&=s\omega+\varepsilon(b_{1}c_{1}\Gamma_{1}(-\alpha)+c_{2}\Gamma_{2}(\theta_{2}-\theta_{1})+b_{2}c_{2}\Gamma_{2}(\theta_{2}-\theta_{3})),\\\
\dot{\theta}_{3}&=s\omega+\varepsilon(b_{1}c_{1}\Gamma_{1}(-\alpha)+c_{2}\Gamma_{2}(\theta_{3}-\theta_{2})),\\\
\dot{\theta}_{4}&=s\omega+\varepsilon(c_{1}\Gamma_{1}(\alpha)+b_{2}c_{2}\Gamma_{2}(\theta_{1}-\theta_{2})),\\\
\dot{\theta}_{5}&=s\omega+\varepsilon(c_{1}\Gamma_{1}(\alpha)+c_{2}\Gamma_{2}(\theta_{2}-\theta_{1})+b_{2}c_{2}\Gamma_{2}(\theta_{2}-\theta_{3})),\\\
\dot{\theta}_{6}&=s\omega+\varepsilon(c_{1}\Gamma_{1}(\alpha)+c_{2}\Gamma_{2}(\theta_{3}-\theta_{2})).\end{aligned}$
(21)

We next consider the symmetry of the phase differences between the ipsilateral
legs, given in Eq. (15). If the initial phases $\\{\theta_{k}(0)\\}_{k=1}^{6}$
satisfy the symmetry (15) and it is maintained during the subsequent
evolution, their phase differences $\varphi_{12}=\theta_{1}-\theta_{2}$,
$\varphi_{32}=\theta_{3}-\theta_{2}$, $\varphi_{45}=\theta_{4}-\theta_{5}$,
and $\varphi_{65}=\theta_{6}-\theta_{5}$ satisfy

$\displaystyle-\varphi_{12}(t)=\varphi_{32}(t)=-\varphi_{45}(t)=\varphi_{65}(t)\coloneqq\beta(t)\mod
2\pi$ (22)

for all $t$. We only need to consider the dynamics of two phase differences
$\varphi_{12}=\theta_{1}-\theta_{2}$ and $\varphi_{32}=\theta_{3}-\theta_{2}$
because $\varphi_{65}=\varphi_{32}$ and $\varphi_{54}=\varphi_{21}$ from Eq.
(17), which are given by

$\displaystyle\dot{\varphi}_{12}$
$\displaystyle=\varepsilon\left(b_{2}c_{2}\Gamma_{2}(\varphi_{12})-c_{2}\Gamma_{2}(-\varphi_{12})-b_{2}c_{2}\Gamma_{2}(-\varphi_{32})\right),$
(23) $\displaystyle\dot{\varphi}_{32}$
$\displaystyle=\varepsilon\left(c_{2}\Gamma_{2}(\varphi_{32})-c_{2}\Gamma_{2}(-\varphi_{12})-b_{2}c_{2}\Gamma_{2}(-\varphi_{32})\right),$
(24)

respectively. Since $-\varphi_{12}=\varphi_{32}=\beta$, the dynamics of the
phase difference $\beta$ can be represented as

$\displaystyle\dot{\beta}=\varepsilon c_{2}\Gamma_{2}(\beta)=-\varepsilon
b_{2}c_{2}\Gamma_{2}(-\beta)$ (25)

from Eqs. (23) and (24).

For $-\varphi_{12}$ and $\varphi_{32}$ to remain equal, Eq. (25) should hold
for all $\beta$, where the PCF $\Gamma_{2}$ should be an odd or even function
because $\Gamma_{2}$ is a continuous nonzero $2\pi$-periodic function and
$c_{2}$ is a non-zero constant. The parameter $b_{2}$ determining excitatory
or inhibitory coupling should be chosen as

$\displaystyle\left\\{\begin{array}[]{ll}b_{2}=1&\quad\mathrm{if\;}\Gamma_{2}\mathrm{\;is\;an\;odd\;function},\\\
b_{2}=-1&\quad\mathrm{if\;}\Gamma_{2}\mathrm{\;is\;an\;even\;function}.\end{array}\right.$
(28)

Thus, the dynamics of the six phases can be described by two independent
equations for the phase differences $\alpha$ between the contralateral legs
and $\beta$ between the adjacent ipsilateral legs,

$\displaystyle\dot{\alpha}$ $\displaystyle=\varepsilon
c_{1}(\Gamma_{1}(\alpha)-b_{1}\Gamma_{1}(-\alpha))\coloneqq\varepsilon\Gamma_{\mathrm{alpha}}(\alpha),$
(29) $\displaystyle\dot{\beta}$ $\displaystyle=\varepsilon
c_{2}\Gamma_{2}(\beta)\coloneqq\varepsilon\Gamma_{\mathrm{beta}}(\beta).$ (30)

By controlling $\alpha$ and $\beta$, we can realize different gait patterns.
The fixed point $(\alpha^{*},\beta^{*})$ of Eqs. (29) and (30) for each gait
is

$\displaystyle(\alpha^{*},\beta^{*})=\left\\{\begin{array}[]{ll}\left(\pi,\frac{1}{3}\pi\right)\quad(\textrm{wave}),\\\
\\\ \left(\frac{2}{3}\pi,\frac{2}{3}\pi\right)\quad(\textrm{tetrapod}),\\\ \\\
\left(\pi,\pi\right)\quad(\textrm{tripod}).\end{array}\right.$ (36)

For these fixed points corresponding to the individual gaits to be realizable,
they should be linearly stable. Moreover, we need the stability of the
periodic solution satisfying the symmetries (14) and (15). These will be
analyzed in the next section.

###  II.4 Stability analysis

We first discuss the stability conditions for the periodic solution of the CPG
network dynamics satisfying the symmetries (14) and (15). We then describe the
conditions for the fixed point $(\alpha^{*},\beta^{*})$ to ensure that the
gaits are stable.

For the periodic solution of the CPG to maintain the symmetry (14), the
variable $\bm{\delta}=\left[\delta_{1}\;\delta_{2}\right]$ defined by

$\displaystyle\delta_{1}$ $\displaystyle=\varphi_{41}-\varphi_{52},$ (37)
$\displaystyle\delta_{2}$ $\displaystyle=\varphi_{63}-\varphi_{52}$ (38)

should be stable at the fixed point $\bm{\delta}=[0\;0]\mod 2\pi$, where
$\varphi_{41}=\theta_{4}-\theta_{1}$, $\varphi_{52}=\theta_{5}-\theta_{2}$,
and $\varphi_{63}=\theta_{6}-\theta_{3}$. Hereafter in this section, for
simplicity, we omit writing ‘mod $2\pi$’ for the fixed points. From Eq. (18),
the dynamics of $\delta_{1,2}$ are represented as

$\displaystyle\begin{aligned}
\dot{\delta}_{1}={}&\dot{\varphi}_{41}-\dot{\varphi}_{52}\\\
={}&c_{1}\Gamma_{1}(\varphi_{41})-b_{1}c_{1}\Gamma_{1}(-\varphi_{41})\\\
&-c_{1}\Gamma_{1}(\varphi_{52})+b_{1}c_{1}\Gamma_{1}(-\varphi_{52}),\end{aligned}$
(39)

and

$\displaystyle\begin{aligned}
\dot{\delta}_{2}={}&\dot{\varphi}_{63}-\dot{\varphi}_{52}\\\
={}&c_{1}\Gamma_{1}(\varphi_{63})-b_{1}c_{1}\Gamma_{1}(-\varphi_{63})\\\
&-c_{1}\Gamma_{1}(\varphi_{52})+b_{1}c_{1}\Gamma_{1}(-\varphi_{52}),\end{aligned}$
(40)

respectively. The fixed point $\bm{\delta}=0$ is linearly stable if all
eigenvalues of the Jacobian $\partial\dot{\bm{\delta}}/\partial\bm{\delta}$
have negative real parts at $\bm{\delta}=0$. By using the chain rule of
derivatives, we obtain

$\displaystyle\begin{aligned}
\frac{\partial\dot{\bm{\delta}}}{\partial\bm{\delta}}&=\begin{bmatrix}\frac{\partial\dot{\delta}_{1}}{\partial\delta_{1}}&\frac{\partial\dot{\delta}_{1}}{\partial\delta_{2}}\\\
\frac{\partial\dot{\delta}_{2}}{\partial\delta_{1}}&\frac{\partial\dot{\delta}_{2}}{\partial\delta_{2}}\end{bmatrix}\\\
&=\begin{bmatrix}\frac{\partial\dot{\delta}_{1}}{\partial\varphi_{41}}-\frac{\partial\dot{\delta}_{1}}{\partial\varphi_{52}}&\frac{\partial\dot{\delta}_{1}}{\partial\varphi_{63}}-\frac{\partial\dot{\delta}_{1}}{\partial\varphi_{52}}\\\
\frac{\partial\dot{\delta}_{2}}{\partial\varphi_{41}}-\frac{\partial\dot{\delta}_{2}}{\partial\varphi_{52}}&\frac{\partial\dot{\delta}_{2}}{\partial\varphi_{63}}-\frac{\partial\dot{\delta}_{2}}{\partial\varphi_{52}}\end{bmatrix}\\\
&=\varepsilon(c_{1}\Gamma^{\prime}_{1}(\alpha)+b_{1}c_{1}\Gamma^{\prime}_{1}(-\alpha))\begin{bmatrix}2&1\\\
1&2\end{bmatrix}.\end{aligned}$ (41)

Therefore, the eigenvalues of $\partial\dot{\bm{\delta}}/\partial\bm{\delta}$
are negative if

$\displaystyle\varepsilon(c_{1}\Gamma^{\prime}_{1}(\alpha)+b_{1}c_{1}\Gamma^{\prime}_{1}(-\alpha))<0.$
(42)

If we use an odd or even function as $\Gamma_{1}$, this condition can be
rewritten as

$\displaystyle\varepsilon c_{1}\Gamma^{\prime}_{1}(\alpha)<0$ (43)

by taking $b_{1}$ as

$\displaystyle\left\\{\begin{array}[]{ll}b_{1}=1&\quad\mathrm{if\;}\Gamma_{1}\mathrm{\;is\;an\;odd\;function},\\\
b_{1}=-1&\quad\mathrm{if\;}\Gamma_{1}\mathrm{\;is\;an\;even\;function}.\end{array}\right.$
(46)

Since $\varepsilon>0$ and $c_{1}>0$, $\bm{\delta}=\bm{0}$ is stable if

$\displaystyle\Gamma^{\prime}_{1}(\alpha)<0.$ (47)

This condition cannot be satisfied for the whole range of $\alpha\in[0,2\pi)$,
but it is sufficient if this condition holds in an interval
$\alpha\in[\alpha_{\mathrm{init}},\alpha^{*}]$ (or
$[\alpha^{*},\alpha_{\mathrm{init}}]$), where $\alpha_{\mathrm{init}}$ and
$\alpha^{*}$ are the initial value and the fixed point of $\alpha$,
respectively.

From now on, we assume that the above condition is satisfied. For the symmetry
(15), or equivalently (22) of the periodic solution to be kept, the variable
$\delta_{3}$ defined by

$\displaystyle\delta_{3}=\varphi_{12}+\varphi_{32}$ (48)

should be stable at the fixed point $\delta_{3}=0$. By adding Eqs. (23) and
(24) and using Eq. (25), the rate of change $\dot{\delta}_{3}$ at
$\delta_{3}=0$ satisfies

$\displaystyle\begin{aligned}
\dot{\delta}_{3}={}&\varepsilon\left(b_{2}c_{2}\Gamma_{2}(\varphi_{12})-2c_{2}\Gamma_{2}(-\varphi_{12})\right)\\\
&+\varepsilon\left(c_{2}\Gamma_{2}(\varphi_{32})-2b_{2}c_{2}\Gamma_{2}(-\varphi_{32})\right)\\\
={}&\varepsilon c_{2}(-\Gamma_{2}(\beta)-b_{2}\Gamma_{2}(-\beta))\\\
={}&0.\end{aligned}$ (49)

The variable $\delta_{3}=0$ is linearly stable if

$\displaystyle\begin{aligned}
\frac{\partial\dot{\delta}_{3}}{\partial\delta_{3}}={}&\frac{\partial\dot{\delta}_{3}}{\partial\varphi_{12}}+\frac{\partial\dot{\delta}_{3}}{\partial\varphi_{32}}\\\
={}&\varepsilon\left(b_{2}c_{2}\Gamma^{\prime}_{2}(\varphi_{12})+2c_{2}\Gamma^{\prime}_{2}(-\varphi_{12})\right)\\\
&+\varepsilon\left(c_{2}\Gamma^{\prime}_{2}(\varphi_{32})+2b_{2}c_{2}\Gamma^{\prime}_{2}(-\varphi_{32})\right)\\\
={}&3\varepsilon
c_{2}(\Gamma^{\prime}_{2}(\beta)+b_{2}\Gamma^{\prime}_{2}(-\beta))<0,\end{aligned}$
(50)

where we used the chain rule of derivatives. This quantity can be calculated
as

$\displaystyle\frac{\partial\dot{\delta}_{3}}{\partial\delta_{3}}=6\varepsilon
c_{2}\Gamma^{\prime}_{2}(\beta)$ (51)

for both odd and even functions of $\Gamma_{2}$ using the condition (28).
Since $\varepsilon>0$ and $c_{2}>0$, $\delta_{3}=0$ is stable if

$\displaystyle\Gamma^{\prime}_{2}(\beta)<0.$ (52)

Although this condition does not hold for all $\beta\in[0,2\pi)$, it is
sufficient if this condition holds in the interval
$\beta\in[\beta_{\mathrm{init}},\beta^{*}]$ (or
$[\beta^{*},\beta_{\mathrm{init}}]$), where $\beta_{\mathrm{init}}$ and
$\beta^{*}$ are the initial value and the fixed point, respectively.

We have shown the conditions for the linear stability of the symmetric
periodic orbit. We now focus on the stability of the fixed point
$(\alpha^{*},\beta^{*})$. The fixed point $\alpha=\alpha^{*}$ is linearly
stable if the following conditions hold:

$\displaystyle\Gamma_{\mathrm{alpha}}(\alpha^{*})$
$\displaystyle=c_{1}(\Gamma_{1}(\alpha^{*})-b_{1}\Gamma_{1}(-\alpha^{*}))=0,$
(53) $\displaystyle\Gamma^{\prime}_{\mathrm{alpha}}(\alpha^{*})$
$\displaystyle=c_{1}(\Gamma^{\prime}_{1}(\alpha^{*})+b_{1}\Gamma^{\prime}_{1}(-\alpha^{*}))<0.$
(54)

If we use an odd or even function for $\Gamma_{1}$ and choose $b_{1}$ as given
by the condition (46),

$\displaystyle\Gamma_{1}(\alpha^{*})$ $\displaystyle=0,$ (55)
$\displaystyle\Gamma^{\prime}_{1}(\alpha^{*})$ $\displaystyle<0$ (56)

are the conditions for the linear stability. The stability condition (56) for
the fixed point $\alpha^{*}$ is consistent with the condition (47) for the
symmetry (14).

The fixed point $\beta=\beta^{*}$ is linearly stable if the following
conditions hold:

$\displaystyle\Gamma_{2}(\beta^{*})$ $\displaystyle=0,$ (57)
$\displaystyle\Gamma^{\prime}_{2}(\beta^{*})$ $\displaystyle<0,$ (58)

which follow from $c_{2}>0$ and Eq. (30). The stability condition (58) for the
fixed point $\beta^{*}$ is consistent with the condition (52) for the symmetry
(15).

##  III Gait control methods

In this section, we explain the methods for gait control. Each gait is
generated by moving the legs when the corresponding phase values exceed a
threshold, where the timing is adjusted by a mutually synchronized solution of
the CPGs with appropriate PCFs $\Gamma_{1,2}$. Gait transitions can be
conducted by changing the PCFs $\Gamma_{1,2}$ between odd and even functions
to switch the fixed points and by alternating $b_{1,2}$ according to the
conditions (46) and (28). The timescale parameter $s$ is adjusted so that the
duration of the swing phase is the same for all of the wave, tetrapod, and
tripod gaits.

###  III.1 Determining gaits from CPG outputs

We first describe the method for determining the swing or stance phase from
the CPG output. We choose one component $x_{\mathrm{out}}$ of the CPG unit
(oscillator) as the output, although the oscillator itself is
multidimensional. The corresponding leg is in the swing phase if the output
$x_{\mathrm{out}}$ is above a threshold $\sigma$ and in the stance phase
otherwise Yu2016gait , i.e.,

$\displaystyle\left\\{\begin{array}[]{ll}\mathrm{swing}&\mathrm{if}\quad
x_{\mathrm{out}}>\sigma,\\\
\mathrm{stance}&\mathrm{otherwise}.\end{array}\right.$ (61)

We assume the output $x_{\mathrm{out}}$ to have only two extremum points
within one period. The threshold $\sigma$ for each gait is determined to
realize the following duty factor $r_{\mathrm{DF}}$:

$\displaystyle
r_{\mathrm{DF}}=\frac{T_{\mathrm{st}}}{T_{\mathrm{sw}}+T_{\mathrm{st}}},$ (62)

where $T_{\mathrm{sw}}$ and $T_{\mathrm{st}}$ represent the durations of the
swing phase and stance phase in one gait period $T/s$, respectively, where
$T_{\mathrm{sw}}+T_{\mathrm{st}}=T/s$ holds. The duty factors for the wave,
tetrapod, and tripod gaits are chosen as $5/6$, $2/3$, and $1/2$,
respectively.

The timescale parameter $s$ is chosen as

$\displaystyle
s=\left\\{\begin{array}[]{ll}\frac{1}{6}\quad(\textrm{wave}),\\\ \\\
\frac{1}{3}\quad(\textrm{tetrapod}),\\\ \\\
\frac{1}{2}\quad(\textrm{tripod}),\end{array}\right.$ (68)

so that the duration of the swing phase of each leg is the same for all of the
wave, tetrapod, and tripod gaits.

###  III.2 Functional forms of the coupling functions

We here describe the odd function $\Gamma_{\mathrm{odd}}$ and even function
$\Gamma_{\mathrm{even}}$ that we use for gait transitions. Both of these
functions are necessary, because for the tripod gait with the fixed point
$(\alpha^{*},\beta^{*})=(\pi,\pi)$, any $2\pi$-periodic odd function
$\Gamma_{\mathrm{odd}}(\varphi)$ has a fixed point at $\varphi=\pi$,
preventing the tripod gait from switching to other gaits, and any
$2\pi$-periodic even function $\Gamma_{\mathrm{even}}(\varphi)$ satisfies
$\Gamma^{\prime}_{\mathrm{even}}(\pi)=0$, where the linear stability of
$\varphi=\pi$ is not guaranteed, i.e., tripod gait is not linearly stable. We
thus switch the odd and even functions to realize gait transitions.

Since the fixed points of the three gaits are $(\alpha^{*},\beta^{*})=$
$(\pi,\pi/3)$, $(2\pi/3,2\pi/3)$, and $(\pi,\pi)$, there can be six cases of
transitions:

$\displaystyle\begin{matrix}\frac{1}{3}\pi\to\frac{2}{3}\pi,&\frac{2}{3}\pi\to\frac{1}{3}\pi,\\\
\\\ \frac{1}{3}\pi\to\pi,&\pi\to\frac{1}{3}\pi,\\\ \\\
\frac{2}{3}\pi\to\pi,&\pi\to\frac{2}{3}\pi.\end{matrix}$ (69)

Here, we choose the odd function as the PCFs $\Gamma_{1,2}$ for the
transitions $\pi/3\to\pi$ and $2\pi/3\to\pi$ because the $2\pi$-periodic odd
function always has a fixed point at $\pi$. We choose the even function for
other transitions.

We next present the odd function $\Gamma_{\mathrm{odd}}$ and even function
$\Gamma_{\mathrm{even}}$ used in this study. We use an odd function defined by

$\displaystyle\Gamma_{\mathrm{odd}}(\varphi)=10\sum_{k=1}^{10}k\exp\left(-\frac{k^{2}}{2}\right)\sin(k\varphi),$
(70)

which is shown in Fig. 3(a). This odd function has only one stable fixed point
at $\varphi^{*}=\pi$ in $(-\pi,\pi]$, which satisfies the stability condition
(56). Also, this odd function satisfies
$\Gamma^{\prime}_{\mathrm{odd}}(\varphi)<0$ in $\varphi\in[\pi/3,\pi]$ for the
stability conditions (47) and (52). We define the even function by

$\displaystyle\Gamma_{\mathrm{even}}(\varphi)=2\cos(\varphi)+1,$ (71)

which is shown in Fig. 3(b). This even function has only one stable fixed
point at $\varphi^{*}=2\pi/3$ in $(-\pi,\pi]$, which satisfies the stability
condition (58). Also, this even function satisfies
$\Gamma^{\prime}_{\mathrm{even}}(\varphi)<0$ in $\varphi\in[\pi/3,\pi]$ for
the stability conditions (47) and (52). We note that, if we consider the
transitions between other gaits as given in Appendix A, we need to design
other even functions depending on the target fixed point.

Figure 3:  (a) Odd function $\Gamma_{\mathrm{odd}}$. (b) Even function
$\Gamma_{\mathrm{even}}$. In each figure, the black dot shows the stable fixed
point and black circle shows the unstable fixed point.

###  III.3 Gait generation and gait transitions

We here explain the method for determining the PCF $\Gamma_{1,2}$ and
$b_{1,2}$ for the gait transitions of the wave, tetrapod, and tripod gaits.

####  III.3.1 Wave gait

Since the fixed point is $(\alpha^{*},\beta^{*})=(\pi,\pi/3)$, we let
$\Gamma_{1}$ be the odd function and $\Gamma_{2}$ be the even function:

$\displaystyle\begin{aligned}
\Gamma_{1}(\alpha)&=\Gamma_{\mathrm{odd}}(\alpha),\\\
\Gamma_{2}(\beta)&=-\Gamma_{\mathrm{even}}(\beta+\pi),\end{aligned}$ (72)

where $\Gamma_{2}(\beta)$ has a stable fixed point at $\beta^{*}=\pi/3$. The
associated parameters $b_{1,2}$ are determined by

$\displaystyle\begin{aligned} b_{1}&=1,\\\ b_{2}&=-1\end{aligned}$ (73)

from the conditions (46) and (28).

####  III.3.2 Tetrapod gait

Since the fixed point is $(\alpha^{*},\beta^{*})=(2\pi/3,2\pi/3)$ in this
case, we choose both $\Gamma_{1,2}$ as the even functions:

$\displaystyle\begin{aligned}
\Gamma_{1}(\alpha)&=\Gamma_{\mathrm{even}}(\alpha),\\\
\Gamma_{2}(\beta)&=\Gamma_{\mathrm{even}}(\beta).\end{aligned}$ (74)

The associated parameters $b_{1,2}$ are determined by

$\displaystyle\begin{aligned} b_{1}&=-1,\\\ b_{2}&=-1\end{aligned}$ (75)

from the conditions (46) and (28).

####  III.3.3 Tripod gait

Since the fixed point is $(\alpha^{*},\beta^{*})=(\pi,\pi)$, we choose both
$\Gamma_{1,2}$ as the odd functions:

$\displaystyle\begin{aligned}
\Gamma_{1}(\alpha)&=\Gamma_{\mathrm{odd}}(\alpha),\\\
\Gamma_{2}(\beta)&=\Gamma_{\mathrm{odd}}(\beta).\end{aligned}$ (76)

The associated parameters $b_{1,2}$ are determined by

$\displaystyle\begin{aligned} b_{1}&=1,\\\ b_{2}&=1\end{aligned}$ (77)

from the conditions (46) and (28).

##  IV Results of Gait Transitions

In this section, we show the results of gait transitions in our CPG network by
numerical simulations.

###  IV.1 FitzHugh-Nagumo CPG model

In our CPG network, various types of limit cycle oscillators can be used as
the CPG unit. We employ the FitzHugh-Nagumo (FHN) oscillator
Fitzhugh1961impulses ; Nagumo1962Active  as the CPG unit, which is a simple
mathematical model of a spiking neuron. The FHN oscillator is described by a
two-dimensional vector $\bm{X}=[x\;y]^{\top}$, which obeys the vector field

$\displaystyle\bm{F}(\bm{X})=d\begin{bmatrix}x-ax^{3}-y\\\
(x+b)c\end{bmatrix},$ (78)

where we assume $(a,b,c,d)=(1/3,0.25,0.15,40)$. The natural period and
frequency of this FHN oscillator without scaling ($s=1$) are $T=0.548$ and
$\omega=11.5$, respectively.

We show the limit-cycle trajectory in Fig. 4(a), each component of the limit-
cycle solution in Fig. 4(b), and the PSF calculated by the adjoint equation in
Fig. 4(c). We choose the first component as the output of the CPG unit, i.e.,
$x_{\mathrm{out}}=x$. The threshold values $\sigma$ for the wave, tetrapod,
and tripod gaits are chosen as $1.6485$, $0.9437$, and $-0.7402$,
respectively.

Figure 4:  Limit cycle and PSF of the FHN oscillator. (a) Limit cycle in the
$xy$ plane. (b) Limit-cycle solution $\bm{X}_{0}=[X_{0,x}\;X_{0,y}]^{\top}$
vs. phase $\theta$. (c) PSF $\bm{Z}=[Z_{x}\;Z_{y}]^{\top}$ vs. phase $\theta$.

###  IV.2 Gait transitions

####  IV.2.1 Wave-tetrapod-tripod

We first simulate the gait transitions from wave to tetrapod and tetrapod to
tripod in order, where the locomotion speed increases. The initial gait is the
wave gait and it is switched to the tetrapod gait at $t=12T_{\mathrm{sw}}$, at
which the PCFs $\Gamma_{1,2}$ and the binary parameters $b_{1,2}$ are changed.
Then, the tetrapod gait is switched to the tripod gait at
$t=36T_{\mathrm{sw}}$, at which $\Gamma_{1,2}$ and $b_{1,2}$ are changed. The
oscillation timescale $s$ is adjusted according to the rule (68). The PCFs
$\Gamma_{1,2}$ and $b_{1,2}$ are changed as described in Sec. III C.

We performed numerical simulations of the dynamics for the phase differences,
Eqs. (29), (30), and the original CPG network dynamics, Eq. (1), for the gait
transitions by using the fourth-order Runge-Kutta method with the time
interval $\Delta_{t}=1.0\times 10^{3}$, where the coupling parameters are
chosen as $c_{1}=4$, $c_{2}=8$, and $\varepsilon=0.1$. We show the outputs of
the CPG units in Fig. 5(a). We can observe that the amplitudes of the outputs
are maintained during the gait transitions, which does not hinder the swing-
stance determination by the threshold value.

The gait diagrams are shown in Fig. 5(b). The purple bars show the wave gait,
where the legs (RH), (RM), (RF), (LH), (LM), and (LF) are moved in order with
the phase shift $\pi/3$. The gray bars between the purple and cyan bars show
the irregular gait during the transition from the wave gait to the tetrapod
gait. The cyan bars show the tetrapod gait, where the three pairs of legs (LH,
RM), (LM, RF), and (LF, RH) are moved in order with the phase shift $2\pi/3$.
The transition time from the wave to tetrapod gait is about
$T_{\mathrm{tr}_{1}}=6.62\times T_{\mathrm{sw}}$. The gray bars between the
cyan and green bars show the irregular gait during the transition from the
tetrapod to tripod gait. Finally, the green bars show the tripod gait, where
the two triplets of legs (RF, LM, RH) and (LF, RM, LH) are moved alternatively
with the phase shift $\pi$. The transition time from the tetrapod gait to the
tripod gait is about $T_{\mathrm{tr}_{2}}=6.13\times T_{\mathrm{sw}}$.

We show the dynamics of the phase differences $\alpha$ and $\beta$ in Figs.
5(c) and (d), respectively. The dynamics of $\alpha$ obtained from Eq. (29)
and the phase difference $\theta_{k+3}-\theta_{k}\;(k=1,2,3)$ evaluated from
the oscillator states by simulating Eq. (1) are slightly different during the
gait transition from $\alpha=\pi$ to $\alpha=2\pi/3$, but both of them
eventually converge to the same target value. There is also a small difference
between $\beta$ obtained from Eq. (30) and the phase differences
$\theta_{2}-\theta_{1}$, $\theta_{3}-\theta_{2}$, $\theta_{5}-\theta_{4}$, and
$\theta_{6}-\theta_{5}$ obtained from the oscillator states obeying Eq. (1).
This is because the reduced phase dynamics and original oscillator dynamics
subjected to weak perturbations are slightly different, causing a small
breakdown of the symmetries, but the phases eventually converge to the target
phase and recover the symmetries because of the stability. Thus, the present
CPG network can achieve the wave-tetrapod-tripod gait transitions by using the
approximate phase-difference dynamics as designed.

Figure 5:  Results for wave-tetrapod-tripod gait transitions. (a) Output
signals of the CPG network. The six signals correspond to the individual legs
as shown in the legend. (b) Gait diagrams over time determined by the rule
(61). The white area shows the stance phase and the colored bars show the
swing phase. The purple bars show the wave gait, the cyan bars show the
tetrapod gait, the green bars show the tripod gait, and the gray bars show the
gaits during transitions, respectively. (c) Dynamics of the phase difference
$\alpha$. The black dashed curve shows the numerically simulated dynamics of
Eq. (29) and the red curves show the dynamics of the phase differences
$\theta_{k+3}-\theta_{k}\;(k=1,2,3)$ measured from the oscillator states
obtained by direct numerical simulations of Eq. (1). (d) Dynamics of the phase
difference $\beta$. The black dashed curve shows the numerically simulated
dynamics of Eq. (30) and the blue curves show the dynamics of the phase
differences $\theta_{2}-\theta_{1}$, $\theta_{3}-\theta_{2}$,
$\theta_{5}-\theta_{4}$, and $\theta_{6}-\theta_{5}$ measured from the
directly-simulated oscillator states obeying Eq. (1).

####  IV.2.2 Tripod-tetrapod-wave

Next, we simulate the tripod to tetrapod and tetrapod to wave transitions in
order, where the locomotion speed decreases. The gait is initially the tripod
gait and it is switched to the tetrapod gait at $t=12T_{\mathrm{sw}}$. Then,
the tetrapod gait is switched to the wave gait at $t=36T_{\mathrm{sw}}$. The
oscillation timescale $s$ is adjusted according to the rule (68). The gait
transition rule is the same as that in the previous case.

We performed numerical simulations of the gait transitions by using the same
procedure as described before. We show the outputs of the CPG network in Fig.
6(a), where the amplitude of the outputs are maintained during the gait
transitions.

We show the gait diagrams in Fig. 6(b). The green bars show the tripod gait,
the cyan bars show the tetrapod gait, and the gray bars between the green and
cyan bars show the irregular gait during the transition from the tripod to
tetrapod gait, respectively. The transition time from the tripod gait to the
tetrapod gait is about $T_{\mathrm{tr}_{1}}=6.62\times T_{\mathrm{sw}}$. From
this and the previous results, we can observe that the tetrapod gait can be
achieved from both the wave and tripod gaits.

The purple bars show the wave gait, and the gray bars between the cyan and
purple bars show the irregular gait during the transition from the tetrapod to
wave gait, respectively. The transition time from the tetrapod to the wave
gait is about $T_{\mathrm{tr}_{2}}=6.13\times T_{\mathrm{sw}}$.

We show the evolution of the phase differences $\alpha$ and $\beta$ in Figs.
6(c) and (d), respectively. The dynamics of $\alpha$ obtained from Eq. (29)
and the phase difference $\theta_{k+3}-\theta_{k}\;(k=1,2,3)$ evaluated from
the oscillator states by simulating Eq. (1) are slightly different during the
gait transitions, but they eventually converge to the same target value. There
are also slight differences between the dynamics of $\beta$ obtained from Eq.
(30) and the phase differences $\theta_{2}-\theta_{1}$,
$\theta_{3}-\theta_{2}$, $\theta_{5}-\theta_{4}$, and $\theta_{6}-\theta_{5}$
evaluated from the oscillator states by simulating Eq. (1), but all of them
eventually converge to the same target value. This is due to the slight
difference between the reduced phase dynamics and original oscillator dynamics
subjected to weak perturbations, which causes a small breakdown of the
symmetries, but the phases eventually go to the target phase and recover the
symmetries from the stability. Thus, the present CPG network can also achieve
the tripod-tetrapod-wave gait transitions.

Figure 6:  Results for tripod-tetrapod-wave gait transitions. (a) Output
signals of the CPG network. The six signals correspond to the individual legs
as shown in the legend. (b) Gait diagrams over time determined by the rule
(61). The white area shows the stance phase and the colored bars show the
swing phase. The green bars show the tripod gait, the cyan bars show the
tetrapod gait, the purple bars show the wave gait, and the gray bars show the
gaits during transitions. (c) Dynamics of the phase difference $\alpha$. The
black dashed curve shows the numerically simulated dynamics of Eq. (29) and
the red curves show the dynamics of the phase differences
$\theta_{k+3}-\theta_{k}\;(k=1,2,3)$ measured from the oscillator states
obtained by direct numerical simulations of Eq. (1), respectively. (d)
Dynamics of the phase difference $\beta$. The black dashed curve shows the
numerically simulated dynamics of Eq. (30) and the blue curves show the
dynamics of the phase differences $\theta_{2}-\theta_{1}$,
$\theta_{3}-\theta_{2}$, $\theta_{5}-\theta_{4}$, and $\theta_{6}-\theta_{5}$
measured from the oscillator states obtained by direct numerical simulations
of Eq. (1).

##  V Conclusions

In this study, we presented a simple CPG network that can control the gait
transitions in hexapod between the different types of gaits whose symmetries
are described in Sec. II B. In our CPG network, the dynamics of the six CPG
oscillators can be reduced to the dynamics of two independent phase
differences between the contralateral legs on the left and right sides and
between the adjacent ipsilateral legs on both sides. By appropriately assuming
the functional form of the mutual coupling functions, the desired
synchronization dynamics can be achieved irrespective of the detailed
properties of the oscillators. Also, since the coupling between the
oscillators is assumed weak, the output signals of the CPGs maintain
approximately constant amplitudes and do not hinder the determination of the
swing-stance phase. We demonstrated that our CPG network can achieve the wave-
tetrapod-tripod gait transitions (speed increasing) and tripod-tetrapod-wave
gait transitions (speed decreasing) smoothly by numerical simulations. The
present CPG network is also applicable to other gait patterns in Appendix A,
which possess the symmetries described in Sec. II B.

The advantage of the present CPG framework is the ability to provide a general
CPG network solely using two independent phase-difference dynamics. Therefore,
the performance of the gait transitions can easily be optimized, compared to
most existing studies that rely on specific models of the CPG oscillators. At
the same time, since the present framework allows us to use various types of
limit-cycle oscillators as the CPG unit, we can employ not only existing
models of oscillators but also artificially designed oscillators
Isjpeert2013dynamical ; Ajalooeian2013general ; Grzelczyk2016prototype ;
Namura2023designing ; Namura2023design  as the CPG unit. In our future work,
we will focus on optimizing the coupling functions within the present CPG
network to achieve faster gait transitions.

###### Acknowledgements.

We acknowledge financial support from JSPS KAKENHI (Nos. JP22K11919 and
JP22H00516) and JST CREST (No. JPMJCR1913).

##  Appendix A Gait variations

In the main text, we considered three typical types of symmetric gaits, i.e.,
the wave, tetrapod, and tripod gaits. However, there exist other gait patterns
that are also called the “wave” or “tetrapod” gait in the literature. We here
list them up to the best of our knowledge. Also, we describe other gait
patterns than the wave, tetrapod, and tripod gaits. We show that all gaits in
this Appendix possess the symmetries described in Sec. II B, which can be
controlled by using the two independent phase differences of our CPG network.

###  A.1 Wave gaits

There are many gait patterns called the “wave” gait, such as (i) (LF), (RF),
(LM), (RM), (LH), (RH) Yu2016gait , (ii) (LF), (LM), (LH), (RF), (RM), (RH)
Bal2021neural  (also called the metachronal tripod in Ref. Inagaki2006wave  or
metachronal+ in Ref. Golubitsky1998modular ), and (iii) (LF), (RM), (LH),
(RF), (LM), (RH) Minati2018versatile  (also called the reversed rolling tripod
in Ref. Inagaki2006wave  or rolling tripod- in Ref. Golubitsky1998modular ),
where (iv) (LF), (RH), (LM), (RF), (LH), (RM) is considered as a similar gait
to (iii), which is called the rolling tripod in Ref. Inagaki2006wave  or
rolling tripod+ in Ref. Golubitsky1998modular . We show the phase shifts when
the phase of LF is $\theta_{1}=0$ and fixed points $(\alpha^{*},\beta^{*})$
corresponding to these gaits in Table 1.

Table 1: Other patterns of wave gaits. | gait patterns | $(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6})$ | $(\alpha^{*},\beta^{*})$
---|---|---|---
(i): | (LF), (RF), (LM), (RM), (LH), (RH) | $\left(0,\frac{4}{3}\pi,\frac{2}{3}\pi,\frac{5}{3}\pi,\pi,\frac{1}{3}\pi\right)$ | $\left(\frac{5}{3}\pi,\frac{4}{3}\pi\right)$
(ii): | (LF), (LM), (LH), (RF), (RM), (RH) | $\left(0,\frac{5}{3}\pi,\frac{4}{3}\pi,\pi,\frac{2}{3}\pi,\frac{1}{3}\pi\right)$ | $\left(\pi,\frac{5}{3}\pi\right)$
(iii): | (LF), (RM), (LH), (RF), (LM), (RH) | $\left(0,\frac{2}{3}\pi,\frac{4}{3}\pi,\pi,\frac{5}{3}\pi,\frac{1}{3}\pi\right)$ | $\left(\pi,\frac{2}{3}\pi\right)$
(iv): | (LF), (RH), (LM), (RF), (LH), (RM) | $\left(0,\frac{4}{3}\pi,\frac{2}{3}\pi,\pi,\frac{1}{3}\pi,\frac{5}{3}\pi\right)$ | $\left(\pi,\frac{4}{3}\pi\right)$

###  A.2 Tetrapod gaits

In the main text, we considered the tetrapod gait given by (LH, RM), (LM, RF),
(LF, RH), but there can be other types of tetrapod gait patterns
Aminzare2018gait , i.e., (i) (LF, RM), (LH, RF), (LM, RH) called the forward
left tetrapod, (ii) (LF, RH), (LM,RF), (LH, RM) called the backward right
tetrapod, and (iii) (LF, RM), (LM, RH), (LH, RF) called the backward left
tetrapod. Moreover, there exists another tetrapod pattern (iv), in which one
or two legs of (LF,RH), (RM), (LH,RF), (LM) Chen2012smooth ;
Minati2018versatile ; Bal2021neural  are moved simultaneously at equal phase
interval $\pi/2$, which corresponds to an intermediate speed gait between the
wave and tetrapod gait considered in the main text Wilson1966insect . We show
the phase shifts when the phase of LF is $\theta_{1}=0$ and the fixed points
$(\alpha^{*},\beta^{*})$ corresponding to these gaits in Table 2.

Table 2: Other patterns of tetrapod gaits. | gait patterns | $(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6})$ | $(\alpha^{*},\beta^{*})$
---|---|---|---
(i): | (LF, RM), (LH, RF), (LM, RH) | $\left(0,\frac{2}{3}\pi,\frac{4}{3}\pi,\frac{4}{3}\pi,0,\frac{2}{3}\pi\right)$ | $\left(\frac{4}{3}\pi,\frac{2}{3}\pi\right)$
(ii): | (LF, RH), (LM, RF), (LH, RM) | $\left(0,\frac{4}{3}\pi,\frac{2}{3}\pi,\frac{4}{3}\pi,\frac{2}{3}\pi,0\right)$ | $\left(\frac{4}{3}\pi,\frac{4}{3}\pi\right)$
(iii): | (LF, RM), (LM, RH), (LH, RF) | $\left(0,\frac{4}{3}\pi,\frac{2}{3}\pi,\frac{2}{3}\pi,0,\frac{4}{3}\pi\right)$ | $\left(\frac{2}{3}\pi,\frac{4}{3}\pi\right)$
(iv): | (LF, RH), (RM), (LH, RF), (LM) | $\left(0,\frac{3}{2}\pi,\pi,\pi,\frac{1}{2}\pi,0\right)$ | $\left(\pi,\frac{3}{2}\pi\right)$

###  A.3 Other gait patterns

Furthermore, we can consider other gait patterns than the wave, tetrapod, and
tripod gaits, e.g., walking with (i) all legs moving at the same time (called
pronk), (ii) with two triplets of legs (LF, LM, LH) and (RF, RM, RH) moved
alternatively (called pace), (iii) with two sets of legs (LF, LH, RF, RH) and
(LM, RM) moved alternatively (called lurch), (iv) with three pairs of legs
(LF, RF), (LM, RM) and (LH, RH) moved in order, where some pause exists
between the pairs (LH, RH) and (LF, RF), which is called inchworm+, (v) with
three pairs of legs (LF, RF), (LH, RH), and (LM, RM) moved in order, where
some pause exists between the pairs (LH, RH) and (LF, RF), which is called
inchworm-, (vi) with three pairs of legs (LF, RF), (LM, RM), and (LH, RH)
moved at equal phase intervals (called caterpillar+), and (vii) with three
pairs of legs (LF, RF), (LH, RH), and (LM, RM) moved at equal phase intervals
(called caterpillar-) Golubitsky1998modular . We show the phase shifts when
the phase of LF is $\theta_{1}=0$ and the fixed points
$(\alpha^{*},\beta^{*})$ corresponding to these gaits in Table 3.

Table 3: Gait patterns other than the wave, tetrapod, and tripod gaits. | gait patterns | $(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6})$ | $(\alpha^{*},\beta^{*})$
---|---|---|---
(i): | (LF,LM,LH,RF,RM,RH) | $\left(0,0,0,0,0,0\right)$ | $\left(0,0\right)$
(ii): | (LF,LM,LH), (RF,RM,RH) | $\left(0,0,0,\pi,\pi,\pi\right)$ | $\left(\pi,0\right)$
(iii): | (LF,LH,RF,RH), (LM,RM) | $\left(0,\pi,0,0,\pi,0\right)$ | $\left(0,\pi\right)$
(iv): | (LF,RF), (LM,RM), (LH,RH) | $\left(0,\frac{5}{3}\pi,\frac{4}{3}\pi,0,\frac{5}{3}\pi,\frac{4}{3}\pi\right)$ | $\left(0,\frac{5}{3}\pi\right)$
(v): | (LF,RF), (LH,RH), (LM,RM) | $\left(0,\frac{1}{3}\pi,\frac{2}{3}\pi,0,\frac{1}{3}\pi,\frac{2}{3}\pi\right)$ | $\left(0,\frac{1}{3}\pi\right)$
(vi): | (LF,RF), (LM,RM), (LH,RH) | $\left(0,\frac{4}{3}\pi,\frac{2}{3}\pi,0,\frac{4}{3}\pi,\frac{2}{3}\pi\right)$ | $\left(0,\frac{4}{3}\pi\right)$
(vii): | (LF,RF), (LH,RH), (LM,RM) | $\left(0,\frac{2}{3}\pi,\frac{4}{3}\pi,0,\frac{2}{3}\pi,\frac{4}{3}\pi\right)$ | $\left(0,\frac{2}{3}\pi\right)$

## References

  * (1) Norman C. Heglund and C. Richard Taylor.  Speed, stride frequency and energy cost per stride: How do they change with body size and gait?  Journal of Experimental Biology, 138(1):301–318, 09 1988.
  * (2) R Strauss and M Heisenberg.  Coordination of legs during straight walking and turning in Drosophila melanogaster.  Journal of comparative physiology. A, Sensory, neural, and behavioral physiology, 167(3), August 1990.
  * (3) D. Graham and H. Cruse.  Coordinated walking of stick insects on a mercury surface.  Journal of Experimental Biology, 92(1):229–241, 06 1981.
  * (4) M. L. Shik, F. V. Severin, and G. N. Orlovsky.  Control of walking and running by means of electrical stimulation of the midbrain.  Biophysics, 11:659–666, 1966.
  * (5) Sten Grillner.  Locomotion in vertebrates: central mechanisms and reflex interaction.  Physiological Reviews, 55(2):247–304, 1975.
  * (6) Eve Marder, Dirk Bucher, David J. Schulz, and Adam L. Taylor.  Invertebrate central pattern generation moves along.  Current Biology, 15(17):R685–R699, Sep 2005.
  * (7) Sten Grillner.  Biological pattern generation: the cellular and computational logic of networks in motion.  Neuron, 52(5):751–766, 2006.
  * (8) Ilya A. Rybak, Natalia A. Shevtsova, Myriam Lafreniere-Roula, and David A. McCrea.  Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion.  The Journal of Physiology, 577(2):617–639, 2006.
  * (9) Shinya Aoi and Kazuo Tsuchiya.  Locomotion control of a biped robot using nonlinear oscillators.  Autonomous Robots, 19(3):219–232, Dec 2005.
  * (10) Dai Owaki and Akio Ishiguro.  A quadruped robot exhibiting spontaneous gait transitions from walking to trotting to galloping.  Scientific Reports, 7(1):277, Mar 2017.
  * (11) Yuichi Ambe, Timo Nachstedt, Poramate Manoonpong, Florentin WÃ¶rgÃ¶tter, Shinya Aoi, and Fumitoshi Matsuno.  Stability analysis of a hexapod robot driven by distributed nonlinear oscillators with a phase modulation mechanism.  In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 5087–5092, 2013.
  * (12) Ludovico Minati, Mattia Frasca, Natsue Yoshimura, and Yasuharu Koike.  Versatile locomotion control of a hexapod robot using a hierarchical network of nonlinear oscillator circuits.  IEEE Access, 6:8042–8065, 2018.
  * (13) Auke Jan Ijspeert, Alessandro Crespi, Dimitri Ryczko, and Jean-Marie Cabelguen.  From swimming to walking with a salamander robot driven by a spinal cord model.  Science, 315(5817):1416–1420, 2007.
  * (14) J. Estremera, J.A. Cobano, and P. Gonzalez de Santos.  Continuous free-crab gaits for hexapod robots on a natural terrain with forbidden zones: An application to humanitarian demining.  Robotics and Autonomous Systems, 58(5):700–711, 2010.
  * (15) Gang Chen, Yang Han, Yuehua Li, Jiatao Shen, Jiajun Tu, Zhicheng Yu, Junrui Zhang, Hao Cheng, Lvyuan Zhu, and Fei Dong.  Autonomous gait switching method and experiments of a hexapod walking robot for Mars environment with multiple terrains.  Intelligent Service Robotics, Feb 2024.
  * (16) Shinkichi Inagaki, Hideo Yuasa, Takanori Suzuki, and Tamio Arai.  Wave CPG model for autonomous decentralized multi-legged robot: Gait generation and walking speed control.  Robotics and Autonomous Systems, 54(2):118–126, 2006.
  * (17) Ricardo Campos, Vitor Matos, and Cristina Santos.  Hexapod locomotion: A nonlinear dynamical systems approach.  In IECON 2010 - 36th Annual Conference on IEEE Industrial Electronics Society, pages 1546–1551, 2010.
  * (18) Weihai Chen, Guanjiao Ren, Jianbin Zhang, and Jianhua Wang.  Smooth transition between different gaits of a hexapod robot via a central pattern generators algorithm.  Journal of Intelligent & Robotic Systems, 67(3):255–270, Sep 2012.
  * (19) Poramate Manoonpong, Ulrich Parlitz, and Florentin WÃ¶rgÃ¶tter.  Neural control and adaptive neural forward models for insect-like, energy-efficient, and adaptable locomotion of walking machines.  Frontiers in Neural Circuits, 7, 2013.
  * (20) Haitao Yu, Haibo Gao, Liang Ding, Mantian Li, Zongquan Deng, and Guangjun Liu.  Gait generation with smooth transition using CPG-based locomotion control for hexapod walking robot.  IEEE Transactions on Industrial Electronics, 63(9):5488–5500, 2016\.
  * (21) Cafer Bal.  Neural coupled central pattern generator based smooth gait transition of a biomimetic hexapod robot.  Neurocomputing, 420:210–226, 2021.
  * (22) Yaguang Zhu, Liang Zhang, and Poramate Manoonpong.  Generic mechanism for waveform regulation and synchronization of oscillators: An application for robot behavior diversity generation.  IEEE Transactions on Cybernetics, 52(6):4495–4507, 2022.
  * (23) Yoshiki Kuramoto.  Chemical Oscillations, Waves, and Turbulence.  Springer, Berlin, 1984.
  * (24) Frank C. Hoppensteadt and Eugene M. Izhikevich.  Weakly Connected Neural Networks.  Springer, New York, 1997.
  * (25) Arthur T. Winfree.  The Geometry of Biological Time.  Springer, New York, 2001.
  * (26) G. Bard Ermentrout and David H. Terman.  Mathematical Foundations of Neuroscience.  Springer, New York, 2010.
  * (27) Hiroya Nakao.  Phase reduction approach to synchronisation of nonlinear oscillators.  Contemporary Physics, 57(2):188–214, Oct. 2016.
  * (28) Takahiro Harada, Hisa-Aki Tanaka, Michael J. Hankins, and István Z. Kiss.  Optimal waveform for the entrainment of a weakly forced oscillator.  Phys. Rev. Lett., 105:088301, Aug 2010.
  * (29) Isuru Dasanayake and Jr-Shin Li.  Optimal design of minimum-power stimuli for phase models of neuron oscillators.  Phys. Rev. E, 83:061916, Jun 2011.
  * (30) Anatoly Zlotnik, Yifei Chen, István Z. Kiss, Hisa-Aki Tanaka, and Jr-Shin Li.  Optimal waveform for fast entrainment of weakly forced nonlinear oscillators.  Phys. Rev. Lett., 111:024102, Jul 2013.
  * (31) Bharat Monga, Dan Wilson, Tim Matchen, and Jeff Moehlis.  Phase reduction and phase-based optimal control for biological systems: a tutorial.  Biological cybernetics, 113(1):11–46, 2019.
  * (32) Yuzuru Kato, Anatoly Zlotnik, Jr-Shin Li, and Hiroya Nakao.  Optimization of periodic input waveforms for global entrainment of weakly forced limit-cycle oscillators.  Nonlinear Dynamics, 105(3):2247–2263, Aug 2021.
  * (33) Shohei Takata, Yuzuru Kato, and Hiroya Nakao.  Fast optimal entrainment of limit-cycle oscillators by strong periodic inputs via phase-amplitude reduction and Floquet theory.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9):093124, Aug. 2021.
  * (34) Sho Shirasaka, Nobuhiro Watanabe, Yoji Kawamura, and Hiroya Nakao.  Optimizing stability of mutual synchronization between a pair of limit-cycle oscillators with weak cross coupling.  Phys. Rev. E, 96:012223, Jul 2017.
  * (35) Nobuhiro Watanabe, Yuzuru Kato, Sho Shirasaka, and Hiroya Nakao.  Optimization of linear and nonlinear interaction schemes for stable synchronization of weakly coupled limit-cycle oscillators.  Phys. Rev. E, 100:042205, Oct 2019.
  * (36) Norihisa Namura and Hiroya Nakao.  Optimal coupling functions for fast and global synchronization of weakly coupled limit-cycle oscillators.  arXiv:2401.05354, 2024.
  * (37) R. M. Ghigliazza and P. Holmes.  A minimal model of a central pattern generator and motoneurons for insect locomotion.  SIAM Journal on Applied Dynamical Systems, 3(4):671–700, 2004.
  * (38) Zahra Aminzare, Vaibhav Srivastava, and Philip Holmes.  Gait transitions in a phase oscillator model of an insect central pattern generator.  SIAM Journal on Applied Dynamical Systems, 17(1):626–671, 2018\.
  * (39) Azamat Yeldesbay, Tibor Tóth, and Silvia Daun.  The role of phase shifts of sensory inputs in walking revealed by means of phase reduction.  Journal of Computational Neuroscience, 44(3):313–339, Jun 2018\.
  * (40) Eric Brown, Jeff Moehlis, and Philip Holmes.  On the phase reduction and response dynamics of neural oscillator populations.  Neural Computation, 16(4):673–715, Apr. 2004.
  * (41) Donald M. Wilson.  Insect walking.  Annual Review of Entomology, 11(Volume 11, 1966):103–122, 1966\.
  * (42) J. J. Collins and Ian Stewart.  Hexapodal gaits and coupled nonlinear oscillator models.  Biological Cybernetics, 68(4):287–298, Feb 1993.
  * (43) Martin Golubitsky, Ian Stewart, Pietro-Luciano Buono, and J.J. Collins.  A modular network for legged locomotion.  Physica D: Nonlinear Phenomena, 115(1):56–72, 1998.
  * (44) Richard FitzHugh.  Impulses and physiological states in theoretical models of nerve membrane.  Biophysical Journal, 1(6):445–466, 1961.
  * (45) J. Nagumo, S. Arimoto, and S. Yoshizawa.  An active pulse transmission line simulating nerve axon.  Proceedings of the IRE, 50(10):2061–2070, 1962.
  * (46) Auke Jan Ijspeert, Jun Nakanishi, Heiko Hoffmann, Peter Pastor, and Stefan Schaal.  Dynamical movement primitives: Learning attractor models for motor behaviors.  Neural Computation, 25(2):328–373, 02 2013.
  * (47) Mostafa Ajallooeian, Jesse van den Kieboom, Albert Mukovskiy, Martin A. Giese, and Auke J. Ijspeert.  A general family of morphed nonlinear phase oscillators with arbitrary limit cycle shape.  Physica D: Nonlinear Phenomena, 263:41–56, 2013.
  * (48) Dariusz Grzelczyk, Bartosz StaÅczyk, and Jan Awrejcewicz.  Prototype, control system architecture and controlling of the hexapod legs with nonlinear stick-slip vibrations.  Mechatronics, 37:63–78, 2016.
  * (49) Norihisa Namura, Tsubasa Ishii, and Hiroya Nakao.  Designing two-dimensional limit-cycle oscillators with prescribed trajectories and phase-response characteristics.  IEEE Transactions on Automatic Control, pages 1–14, 2023.
  * (50) Norihisa Namura and Hiroya Nakao.  Design of limit-cycle oscillators with prescribed trajectories and phase-response properties via phase reduction and floquet theory.  In 2023 62nd IEEE Conference on Decision and Control (CDC), pages 3962–3968, 2023.