<figure><img src="image_1.png" /></figure>

HAL
open science

# Symmetry breaking and gait transition induced by hydrodynamic sensory feedback in an anguilliform swimming robot

Johann Herault, Laura Paez, Kamilo Melo, Robin Thandiackal, Vincent Lebastard, Frédéric Boyer, Auke Ijspeert

► To cite this version:

Johann Herault, Laura Paez, Kamilo Melo, Robin Thandiackal, Vincent Lebastard, et al. Symmetry breaking and gait transition induced by hydrodynamic sensory feedback in an anguilliform swimming robot. Physical Review E, 2024, 110 (5), pp.055104. 10.1103/PhysRevE.110.055104. hal-04817887

HAL Id: hal-04817887

https://hal.science/hal-04817887v1

Submitted on 4 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Symmetry breaking and gait transition induced by hydrodynamic
sensory feedback in an anguilliform swimming robot.

Johann Herault ¹,* Laura Paez ², Kamilo Melo ³, Robin Thandiackal

⁴, Vincent Lebastard ¹, Frédéric Boyer ¹, and Auke Ijspeert ²

¹IMT Atlantique, LS2N, 4 Rue Alfred Kastler, 44300 Nantes, France

²École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

³ KM-RoBoTa, Renens, Switzerland and

⁴Harvard University, Cambridge MA, USA

(Dated: December 4, 2024)
Abstract

The goal of this article is to identify and understand the fundamental role of spatial symmetries
in the emergence of undulatory swimming using an anguilliform robot. Here, the local torque at the
joints of the robot is controlled by a chain of oscillators forming a central pattern generator (CPG).
By implementing a symmetric CPG with respect to the transverse plane, motor activation waves
are inhibited, preventing the emergence of undulatory swimming and resulting in an oscillatory
gait. We show experimentally that the swimmer can recover from the traveling wave inhibition
by using distributed fluid force feedback to modulate the phase dynamics of each oscillator. This
transition from oscillatory to undulating swimming is characterized by a symmetry breaking in the
CPG and the body dynamics. By studying the stability of the oscillator chain, we show that the
sensory feedback produces a frequency detuning gradient along the CPG chain while preserving
its stability. To explain the origin of the instability, we introduce a toy model where the couplings
between the dynamics of the oscillators and the body deformation reinforce the symmetry breaking

# I. INTRODUCTION

Aquatic animals have developed a wide variety of gaits to move efficiently in various fluid environments. Great progress has been made over the last 20 years in understanding these swimming modes by reproducing as closely as possible the interactions between an animal and its environment thanks to bioinspired robots [1–4] or integrative numerical models with realistic hydrodynamic conditions [5–7]. These platforms have shown how complex motion, including gait and maneuvers, results from the interplay among the swimmer morphology, the nervous system and the aquatic environment.

Conversely, the increase in the number of degrees of freedom of these systems makes it
more complex to analyze the fine mechanisms of gait selection, which are generally inter-
preted as emergent phenomena. Nonetheless, these dynamics include also regularities and
spatial symmetries that could help us to reduce this apparent complexity [8, 9]. In the
context of anguiliform swimming, the symmetry of motor activation with respect to the
transversal plane of the swimmer determines the nature of the gait (see Fig. 1). Here, we

<footnote>* johann.heraul@imt-atlantique.fr</footnote>
study experimentally the role of this symmetry to better characterize the interactions be-
tween motor control and swimmer dynamics through sensory feedback, a mechanism which
remains yet poorly understood.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Illustration of the symmetries in swimming lampreys. If the CPG couplings are symmetric with respect to the transversal plane, the resulting wave of body deformation is also symmetric and standing. A necessary condition to drive a propagative wave is to break this symmetry, which results in two patterns of body deformation with distinct symmetries (even and odd).</figcaption></figure>

We will consider undulatory swimming of anguilliform swimmers at large Reynolds num-
bers such as eels, or lampreys. This gait consists in accelerating the surrounding fluid down
the body thanks to a head-to-tail traveling bending wave. In vertebrates, the rhythmic
muscular pattern sustaining this wave is generated by distributed neural oscillators located
in the spinal cord. This complex network of oscillators [10] is called the central pattern
generator (CPG) [11], and is the basic unit of motor control in the central system.

In lampreys, there is clear evidence that the CPG maintains phase lags along the spinal
cord resulting in a traveling wave of muscular activity [12]. Experiments have shown that
a traveling wave of neural activity called *fictive locomotion* [11], can be sustained by only
the spinal cord irrespective of the peripheral nervous system, i.e. without sensory feedback
[13]. Nonetheless, sensory signals provided by the peripheral nervous system play a crucial
role in locomotion by shaping and regulating the rhythmic patterns via stimuli from the
body stretching [14] or the surrounding environment [11, 15]. In salamanders, the open loop
CPG produces all kinds of traveling waves with positive phase and negative phase lags [16]. Hence, the CPGs are not necessarily well-tuned for all animals and sensory feedback could be a key mechanism in the gait selection process in this case.

Sensory feedback, such as fluid response or stretch feedback, can provide better adaptability to the aquatic environment, such as gait selection [17] or optimal beat frequency [18]. Robotic experiments [17] have also demonstrated that undulatory swimming can be generated by pressure force feedback without CPG couplings. These observations support that hydrodynamic sensory feedback coupled with a CPG model reinforces the wave entrainment and increases the robustness against neural disabilities [17]. A recent numerical study [7] has confirmed that a mechanism of compensation can also be obtained by proprioceptive feedback given by the local body curvature. While these results demonstrate the importance of sensory feedback in motor coordination, the explanation and the mathematical modeling of these observations remain unclear.

More specifically, we want to know how sensory feedback contributes to the coordination and recovery of motor function through the prism of symmetries. To study this mechanism, we implement an artificial CPG such that its symmetries prevent the development of traveling neural waves along the oscillator chain. With such a setup, we can now study how sensory feedback can generate a propagating wave in a symmetric network of oscillators. By tackling this new problem, we aim to precisely analyze the mechanism of the gait selection process induced by sensory feedback (see Fig. 2). To do this, we will use the tools of nonlinear physics, in particular synchronization theory, modal approach, and perturbation theory.

This paper is structured as follows. After presenting the robot (section II) and the artificial neural system (section III), we show in section IV that a propagating wave of deformation along the swimmer can spontaneously emerge via a symmetry breaking bifurcation when the strength of the feedback in the network exceeds a threshold. This bifurcation is associated with a gait transition from the (zero-phase lag) oscillatory swimming to the (negative phase lag) undulatory swimming. This observation shows that a feedback loop provided by the sensory system ensures the emergence of an activation wave despite a severe disruption in the CPG [7, 17] (see Fig. 2).

To account for the observed instability, we propose a theoretical framework for quantitatively analyzing the stability of the oscillator network (section V). In this theoretical
<figure><img src="image_4.png" /><figcaption>FIG. 2. Schematic illustration of the role of the symmetries in the process of gait selection. The body dynamics and the fluid response "reflect" the symmetries of motor coordination. Hence, an initial antisymmetric component in the neural system can be reinforced by positive feedback allowed by sensory feedback. Under specific conditions, this loop can be at the origin of an activation wave in a spatially invariant CPG designed to inhibit propagating waves.</figcaption></figure>

framework, we explain how exteroceptive sensory feedback gives rise to a network whose
couplings are determined by the physics of swimming. Here, the physical medium for com-
munication between the oscillators is the momentum exchanged between the servomotors
and the sensors. We then explain the shape of this network, showing that these couplings
are linked to the transfer function between torque and sensory feedback. We show how to
analyze the stability of this network using a perturbative approach. The main result of this
theoretical part is that sensory feedback does not change the stability of the oscillator net-
work and produces a frequency detuning gradient along the oscillator chain. This approach
gives a good quantitative agreement with the experimental results. Finally, we propose a
simple toy model to account for the instability based on reciprocal couplings between the
oscillators and the body dynamics.
## II. THE ROBOTIC PLATFORM

### A. AgnathaX

We investigate the role of symmetry in the CPG thanks to a modular lamprey like robot, called AgnathaX (see Fig. 3), built out of ten active body modules, a head unit and a passive tail (more details can be found in [17]). Each of the ten body modules is actuated by a servomotor. The onboard computer acquires the external sensor signals from all body modules as well as the full state of the motors (i.e. position, speed, voltage, and current). It also evaluates the distributed controllers for each module in a 100Hz control loop. The corresponding output commands of each controller are sent back to the motors to move the robot. Data gathered from the sensors, motors, and the states of the controllers are logged during swimming experiments and is exploited for posterior analysis.

The shape of the tail is designed to ensure that it resonates at nominal robot undulation angles of ±30°, and a frequency of 0.75 Hz [19]. A protective soft, hydrophobic and highly flexible wrapping sleeve covered the robot modules.

A feed-forward controller based on current control is designed [17] to control the torque exerted on the *i*th servomotor. This torque is then proportional to an activation signal $a_i = \cos(\phi_i)$ with a time-varying phase $\phi_i(t)$ (see Fig. 3) resulting from the dynamics of the *i*th oscillator in the oscillator chain whose model is introduced in the next section (Sec. III).

<figure><img src="image_5.png" /><figcaption>FIG. 3. (color online) The robot without (a) and with its suit (b). (c) Local differential hydrodynamic force sensors. (d) Schema of the phase dynamics where each oscillator is coupled to nearest oscillators and modulated by the local hydrodynamic force.</figcaption></figure>
## B. Sensory feedback

The robot is equipped with load cells reproducing the mechanoreceptors on the animal skin. In lampreys, they are called dorsal cells, and are sensitive to touch and pressure on the lateral side of the body [20]. Hence, pairs of contralateral sensing modules located on both sides of each segment measure the differential lateral hydrodynamic forces $F_i(t)$ exerted on each body section [17] via two load cells. For convenience, we prefer to use the nondimensionalized form $S_i(t) = F_i(t)/\tilde{F}$ with $\tilde{F} = 1N$ to simplify the principle of homogeneity of the equations. These data provide the exteroceptive feedback, which accounts for the response of the environment to the body motion.

# III. THE CPG MODEL

## A. Phase equation

To drive the motor, we use a motor control based on central pattern generators (CPG). Many models can reproduce the main features of the CPG such as dynamical systems coupled to each other or multiple arrays of oscillators. Here, we consider a simple model of a 1D array of oscillators such that the dynamics of the ith oscillator is governed by its phase $\phi_i(t)$. The local oscillator drives locally the bending torque that is normally produced by antagonist muscles, each one being activated by one neural oscillator. By using this 1D model, we focus only on the mechanism of synchronization along the rostrocaudal axis as well as the symmetries with respect to the transverse plane.

For the sake of simplicity, we also assume that the chain features only close neighbor couplings as suggested by previous studies [2, 21], and each phase is modulated by local sensory feedback $S_i(t)$. Based on a standard phase reduction [22], the time evolution of the phase $\phi_i(t)$ can be written under the very general form

$$ \dot{\phi}_i = \omega_i + \sum_{k=\pm 1} C_{i,i+k} \sin(\phi_{i+k} - \phi_i + \psi_{i,i+k}) + Z_i(\phi_i)S_i(t), \quad (1) $$

where $i \in [1, N]$, N is the index of each oscillator, $\omega_i$ the intrinsic frequency of the local oscillator, $C_{i,k}$ and $\psi_{i,i+k}$ the coupling strength and phase delay, and finally $Z_i(\phi_i)$ the sensitivity function that accounts for the response of the local oscillator to exogenous perturbations.
The second right-hand side (RHS) term represents the interoscillator couplings. The last RHS term corresponds to the contribution of the peripheral system to the synchronization process, while the other terms are related to the CPG configuration.

Various strategies lead to propagating waves by tuning the CPG properties [9, 11, 21]. In the absence of sensory feedback, the main mechanisms are (i) a phase delay $\psi_{i,i+k}$ along the spinal cord, (ii) an inhomogeneous distribution of intrinsic frequency $\omega_i$ with either a frequency gradient or master-slave oscillators, or (iii) an asymmetry between ascending $C_{i,i-1}$ and descending $C_{i,i+1}$ couplings. Note that these mechanisms can be combined to enhance the redundancy in the CPG. Conversely, if we consider identical oscillators ($\omega_i = \omega_0$) with symmetrical couplings between close neighbors ($C_{i,i\pm1} = C$) without delay ($\psi_{i,i+k} = 0$), i.e. a spatially invariant system, no phase lag occurs.

Although this latter configuration has never been observed in nature, it represents an interesting paradigm to investigate the effect of sensory feedback on the dynamics of CPG. Indeed, the appearance of the propagating motor activation wave can only be caused by sensory feedback in this case. By pursuing our reasoning on symmetries, this symmetry breaking must also take the form of one of the three proposed mechanisms. Such mechanisms are particularly challenging to identify or characterize if the CPG is not initially symmetric and motivates the implementation of a symmetric CPG.

Consequently, we propose to reduce Eq. (1) to the following simple model

$$ \dot{\phi}_i = \omega_0 + C \sum_{k=\pm 1} \sin(\phi_{i+k} - \phi_i) + \sigma \cos(\phi_i) S_i(t), \quad (2) $$

with a uniform intrinsic frequency $\omega_0 = 0.75 \times 2\pi$ rad/s, and a coupling strength set to $C = 10$ rad/s. The sensitivity function $Z_i(\phi_i)$ is fixed and given by $Z_i(\phi_i) = \sigma \cos(\phi_i)$, which is the same as our previous study [17].

The free parameter $\sigma$ can be varied to modify the weight of sensory feedback in the chain. The differential equation given by Eq. (2) is integrated in real-time in each robot module, such that the torque provided by each servomotor is proportional to the activation signal $\tau_i \propto a_i$ with $a_i = \cos(\phi_i)$. It is worth noting that the motor activation $a_i$ and the sensitivity function are defined by the same trigonometric function $\cos(\phi_i)$. In our previous study [17], we observed that this choice led to the emergence of propagating waves in the network.
## B. Phase decomposition and symmetry

To reduce the complexity of the phase dynamics, we propose to solve the linearized equation (2) without sensory feedback (last RHS term of Eq. (2)) and establish a set of normal modes to reduce the complexity of the phase dynamics. In the limit of a small phase difference $\phi_{i+1} - \phi_i$, equation (2) without sensory feedback can be linearized such that

$$ \dot{\phi}_i = \omega_0 + C(\phi_{i+1} + \phi_{i-1} - 2\phi_i). \qquad (3) $$

This schema is called diffusive coupling [23, 24] and is linearly stable. To show this property, we decompose the phase $\phi_i$ as

$$ \phi_i = \Omega t + \delta\psi_i(t) \qquad (4) $$

with $\Omega$ the mutual frequency of the oscillators and $\delta\psi_i(t)$ the local phase lag, which varies during the setup of the synchronization process, and then converges to a constant value. In the context of motor activation, the angular frequency $\Omega$ controls the beating amplitude. The phase lags $\delta\psi_i$ along the spinal cord determine the nature of the gait, i.e. zero phase lag for the oscillatory swimming [25] and a negative phase lag for the undulatory swimming.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Top: eigenvalues $\hat{C}_n$ as a function of the index $n$. Bottom: eigenmodes with $n = 1$ (blue or black curve), $n = 2$ (red or grey) and $n = 3$ (yellow or light curve).</figcaption></figure>

The diffusive schema is characterized by eigenmodes $\mathbf{v}_n$ with eigenvalues $\hat{C}_n$ defined by

$$ \left\{ \begin{aligned} v_n(i) &= \cos(\kappa_n(i - 0.5)), \\ \hat{C}_n &= -2C(1 - \cos(\kappa_n)), \end{aligned} \right. \qquad (5) $$
with the wave vector $\kappa_n = n\pi/N$ and $n \in [0, N - 1]$. These eigenvalues and eigenmodes (for $n = 1 - 3$) are reported in Fig. 4. The mutual frequency $\Omega$ corresponds to the neutral mode with $n = 0$, associated with the uniform mode $v_0(i) = 1$. Hence, the phase lags $\delta\psi_i(t)$ can be decomposed into

$$ \delta\psi_i(t) = \sum_{n=1}^{N} \delta\hat{\psi}_n(t)v_n(i). \qquad (6) $$

We have introduced the modal component $\delta\hat{\psi}_n(t)$ that characterizes the amplitude of the nth mode. This modal component is computed by modal projection defined as $\delta\hat{\psi}_n = \langle \mathbf{v}_n, \delta\psi \rangle$ and given by

$$ \langle \mathbf{v}_n, \delta\psi \rangle = \frac{1}{\sum_i v_n(i)^2} \sum_i v_n(i) \delta\psi_i. \qquad (7) $$

The diffusive coupling reduces the spatial fluctuation at small scales with a damping rate given by $\hat{C}_n$. By projecting Eq. (3) based on Eq. (7), one obtains an ordinary differential equation for each modal component

$$ \frac{d}{dt}\delta\hat{\psi}_n = \hat{C}_n\delta\hat{\psi}_n. \qquad (8) $$

Consequently, in the absence of sensory feedback, the corresponding modal components $\delta\hat{\psi}_n(t)$ are given by

$$ \delta\hat{\psi}_n(t) = \delta\hat{\psi}_n(0)e^{\hat{C}_n t}. \qquad (9) $$

Similarly, one shows that the mutual frequency is $\omega_0$ such that after a transient, the oscillator phase converges to $\phi_i = \omega_0 t$ with $\delta\psi_i = 0$. We have confirmed that in the absence of feedback, the oscillators are synchronized. The resulting activation wave is standing, as well as the wave of body deformation illustrated in Fig. 5 (top) because all the motors beat in sync. Here, the propagating wave of activation can only arise from sensory feedback terms.

The decomposition proposed in Eq. (4) requires filtering out the fast time fluctuations arising from the nonlinear dynamics, a point treated in Sec. V. In this sense, the synchronization process is only meaningful for the slow dynamics, i.e. on a timescale larger than $2\pi/\Omega$.
<figure><img src="image_7.png" /><figcaption>FIG. 5. Superimposition of pictures illustrating the robot gaits (nearly every seconds). They correspond to the oscillatory swimming (top) and the undulatory swimming (bottom) respectively before and past the threshold of the instability. Oscillatory swimming, for example, can be exploited by thunniform swimmers[25], while slender swimmers such as anguilliform swimmers use undulatory swimming.</figcaption></figure>

The spatial modes {$v_n$} will allow us to reduce the complexity of the network dynamics by projecting the phase dynamics onto a few modes with different spatial symmetries with respect to the transverse plane. The even modes starting with the $n = 2$ mode are symmetric with respect to the transversal plane of the chain, while the odd modes $n \ge 1$ are antisymmetric. The decomposition given by Eq. (4) suggests that the large-scale modes ($n = 1, 2$) are the least damped (cf Fig. 4) and therefore a priori dominant while having distinct symmetries. Hence, our approach focuses on the modal dynamics restricted to a few modes rather than the local dynamics of each oscillator to characterize the appearance of a propagating wave. We will see that these symmetries and the large-scale dynamics are decisive in the gait selection process.

# IV. GAIT TRANSITION INDUCED BY SENSORY FEEDBACK

## A. Gait transition and phase bifurcation

In this section, we report experimental results on the effect of sensory feedback in the oscillator chain on swimming performance. We gradually increase $\sigma$ from 0 to 1 in Eq. (2) with a set of 125 experiments, with at least 5 measurements for each $\sigma$. After a transient, all
the runs converge to a steady regime. For each run, we compute a mean frequency detuning $\delta\Omega = \overline{\dot{\phi}_i - \omega_0}$, with the time averaging operator $\overline{\cdot}$, and a average phase shift $\delta\psi_i = \overline{\phi_i - \Omega t}$ with $\Omega = \omega_0 + \delta\Omega$, and its modal amplitude given by the projection $\delta\psi_n = \langle \mathbf{v}_n, \delta\psi \rangle$. The spatial average of $\delta\psi_i$ is set to zero by shifting the time origin.

<figure><img src="image_8.png" /><figcaption>FIG. 6. Measured (open circle) and predicted (crosses) frequency detuning $\delta\Omega$ as a function of the control parameter $\sigma$. The prediction is based on the solution of the phase equation given by Eq. (18) (section V).</figcaption></figure>

For $\sigma < 0.175$, the oscillators are fully synchronized with $\delta\Omega \simeq 0$ (see Fig. 6) and all the phase shifts $\delta\psi_i$ remain small (dotted lines in Fig. 7 (a)). In this parameter range, the robot performs an oscillatory swimming (Fig. 5, top) with a forward body velocity $U_b \simeq 0.22\,\text{m/s}$ (Fig. 8 (a)). When the parameter $\sigma$ crosses a threshold $\sigma_c = 0.175 \pm 0.05$, the body velocity $U_b$ shows a sharp transition with a significant increase of 85% (Fig. 8 (a)). This is a noticeable performance gain since all the other physical parameters are kept constant. This transition can be attributed to a sudden change of gait as illustrated in Fig.5.

A wave of body deformation can be identified in Figs. 5 (bottom) and 8(c-e). During this transition, the oscillators remain synchronized with a frequency detuning $\delta\Omega$ that increases almost linearly beyond the threshold (Fig. 6) as previously observed [17]. Reciprocally, we observe in the network the progressive establishment of a negative gradient of phase shift $\delta\psi_i$ along the oscillator chain (thick black curves in Fig. 7.(a)) whose amplitude increases with the parameter $\sigma$ (grey arrows). Our results demonstrate that sensory feedback can trigger a sharp transition from oscillatory to undulatory swimming for anguilliform swimmers.

This transition is associated with a symmetry breaking in the network as illustrated by the antisymmetric curve of the phase lag $\delta\psi_i$ (see Fig. 7). Thanks to a modal analysis, we can show that this symmetry breaking comes from an unstable odd mode. Using the
<figure><img src="image_9.png" /><figcaption>FIG. 7. (a) Phase lag $\delta\psi_i$ along the oscillator chain $i \in [1, 10]$ for different $\sigma$. The dashed and dotted curves correspond respectively to $\sigma = 0$ and $0.1$. Above the onset, the thick black curves correspond to $\sigma = [0.18, 0.21, 0.3, 0.5, 0.8, 1]$. (b) Measured (open circle) and predicted (cross) modal components $\delta\hat{\psi}_1$ (blue or dark circles) and $\delta\hat{\psi}_2$ (red or grey diamonds) of the phase lag $\delta\psi_i$ for each $\sigma$ (see Eq. (4)). The prediction is based on the solution of the phase equation given by Eq. (18) (section V).</figcaption></figure>

modal decomposition of the phase shift $\delta\psi_i$, we observe that the amplitude of the first odd component $\delta\hat{\psi}_1 = \langle \mathbf{v}_1, \delta\psi \rangle$ increases significantly after the threshold $\sigma_c$ (Fig. 7 (b)). The trend of this curve is a clear signature of a critical bifurcation with a small imperfection. This odd mode breaks the initial even symmetry of the oscillator chain and the motor coordination. This modal analysis shows that the emergence of a phase gradient in the oscillator chain is explained by the bifurcation of the first eigenmode of the oscillator chain. We have also reported the amplitude of the second eigenmode $\delta\hat{\psi}_2$, which is small enough to be neglected in the phase dynamics.

To further characterize this symmetry breaking, we also computed the growth rate of the instability (cf Fig. 9). The modal phase shift $\delta\hat{\psi}_1$ initially describes an exponential growth before reaching a saturated state. By estimating the slope of the tangent at the origin, one can compute a growth rate $\xi$ ($s^{-1}$), which is then plotted as a function of $\sigma$ for different realizations (black dots). Surprisingly, this value does not increase linearly with $\sigma - \sigma_c$ as expected for a standard pitchfork instability. Instead, we observe that the growth rate has a square root dependence. This observation suggests that the instability is not described by a single mode with one amplitude equation, as in the pitchfork instability, but rather by at least two interacting modes. This point will be important in determining the model that accounts for the instability (section V).
<figure><img src="image_10.png" /><figcaption>FIG. 8. (Color online) (a) Asymptotic forward body velocity $U_b$ of the robot as a function of $\sigma$. (b-e) Superimposition of polynomials interpolation of the body deformation in the barycentric reference frame for different values of $\sigma$. From ochre or light ($t=0$) to dark blue (end of the cycle).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. Top: time evolution of the component $\delta\hat{\psi}_1(t)$ for various $\sigma$ with $\sigma \in [0.18, 0.22, 0.5, 0.9]$ (recognisable by their saturation value which increases with $\sigma$)). The slope of the tangent at the origin determines the initial growth-rate. Bottom: interpolation of the growth-rate $\xi$ as a function of $\sigma$. The red curve corresponds to the critical variation $\xi = -\nu + \sqrt{\nu^2 + \alpha(\sigma - \sigma_c)}$ with only two free parameters $\nu = 0.0165\text{s}^{-1}$ and $\alpha = 0.902\text{s}^{-1}$ (cf. equation (23)).</figcaption></figure>

## B. Characterization of the internal body motion

In this section, we examine the response of the body dynamics during this transition. In figure 10, we show the spatiotemporal dynamics of the joint $\theta_i(t)$ before the instability with $\sigma = 0$ (top), just above the threshold with $\sigma = 0.18$ (middle), and far beyond the
threshold of the instability with $\sigma = 1$ (bottom). Before the instability, the periodic state is almost symmetric with respect to the transverse plane in agreement with the symmetry of the motor activation. When the parameter $\sigma$ is above the threshold $\sigma_c = 0.175$ (Fig. 10, middle), one observes the progressive build up of the wave traveling from the top (the head) to the bottom (the tail) of the image. As $\sigma$ increases, the propagating wave becomes clearer.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Spatiotemporal diagram of the joint angle amplitude (index of the motor versus time) for (a) $\sigma = 0$, (b) $\sigma = 0.18$, and (c) $\sigma = 1$ (c).</figcaption></figure>

According to Curie's principles, the symmetry breaking in the oscillator dynamics must be reflected in the body dynamics via motor activation. It is tempting to assume a correspondence between the modal dynamics in the network and the internal body motion. To test this hypothesis, we compare the most energetic bending modes (thick curves in Fig. 11(a)) obtained by a proper orthogonal decomposition analysis [26] with the first modes $\mathbf{v}_1$ and $\mathbf{v}_2$ (dashed curves in Fig. 11(a)). The eigenmodes of the joint amplitudes look very similar to the first modes of the vector set $\{\mathbf{v}_n\}$ (cf section III B), since they are similar to those of the lumped-mass/stiffness model [27]. We will therefore use this basis to capture the internal dynamics of the body and facilitate the comparison with the modal phase dynamics.

Thus, we introduce the modal decomposition of the joint amplitude components $\theta_i$ given by the modal amplitude $\hat{\theta}_n = \langle \mathbf{v}_n, \boldsymbol{\theta} \rangle$. We compute its mean square amplitude for each $\sigma$ and $n = 0, 1, 2$ in Fig. 11(b). Before the instability, the oscillation is mostly composed of even
modes ($\hat{\theta}_0$ and $\hat{\theta}_2$ in Fig. 11) that are consistent with the symmetry of the driving torques. The small odd component $\hat{\theta}_1$ is mostly due to the asymmetry of the robot (heavier head and softer tail), and could explain the imperfection in the bifurcation. Near the threshold, the amplitudes of the even modes $\hat{\theta}_0$ and $\hat{\theta}_2$ decrease. Conversely, the odd mode $\hat{\theta}_1$ increases significantly (Fig. 11, (b)) which is correlated with the setup of the propagating wave (cf Fig. 8). After the onset, the component $\hat{\theta}_2$ grows again, probably due to the nonlinearities of the motor activation.

<figure><img src="image_13.png" /><figcaption>FIG. 11. (a) Mean square amplitude of the joint modal component $\hat{\theta}_n$ as a function of $\sigma$. (b) Two first modes from POD for $\sigma = 0.6$ (thick blue and red curves) superimposed to the modes $v_1$ and a combination of $v_0$ and $v_2$ (dashed curves).</figcaption></figure>

The emergence of the propagating wave is actually correlated to the presence of modes $\hat{\theta}_1$ and $\hat{\theta}_2$ of different parity. The superposition of an odd and an even mode $\hat{\theta}_1$ and $\hat{\theta}_2$ is a general mechanism that accounts for wave propagation in bounded domain [28, 29]. Furthermore, the bifurcation imperfection caused by the robot asymmetries, which favors an odd mode sign, is also consistent with the reproducible rostrocaudal direction of bending wave propagation for all $\sigma > \sigma_c$ [17].
<figure><img src="image_14.png" /><figcaption>FIG. 12. Analogy between the communication between fireflies made possible by the emission of flashes and the measurement of luminosity, and our experimental platform where each oscillator recovers information on the state of all the oscillators by perceiving the hydrodynamic force acting the corresponding module. The angular momentum injected by the k-th servomotor into the serial mechanical structure is then transferred to the fluid and measured on the wall of the i-th module. The dynamics of the oscillator of this module is therefore modified by the activity of all the servomotors, where the weight of the couplings between oscillators depends on the fluid-structure interactions and the viscoelastic properties of the robot.</figcaption></figure>

# V. STABILITY ANALYSIS

## A. Physical interactions induced a complex network of oscillators

The simultaneous symmetry breaking in the oscillator chain and in the body motion suggests a gait transition driven by a positive feedback between the phases and the joint dynamics. If the coupling from the oscillator chain to the body motion is provided by the activation of the servomotors, the effect of sensory feedback on the phase dynamics is more complex to interpret.

Our first goal is to describe the contribution of sensory feedback $S_i$ to the dynamics of the i-th oscillator. A previous study [30] has highlighted the fact that the gait transition in a CPG can be explained by phase entrainment and phase-locking mechanisms induced by sensory feedback. However, they hypothesized that sensory feedback is independent of phase dynamics $\phi_i$ to develop an analogy with a particle moving in a potential landscape. On the contrary, here we hypothesize that sensory feedback results from the dynamics of all the phases, creating a complex network of oscillators, whose stability must be determined.
To illustrate this hypothesis, we propose to develop an analogy between the phenomenon of synchronization of a population of fireflies, a canonical example of collective synchronization, and our problem. This comparison is illustrated in figure 12. Male fireflies synchronize their emission of light of intensity $I_i$ in order to optimize their chances of meeting a mate [31]. The intensity modulation is controlled by the phase of an oscillator $\phi_i$, which can be written as $I_i \propto \cos(\phi_i)$ for the purpose of the analogy. This synchronization is only possible if each firefly perceives the local luminosity $S_i$ resulting from the flashing of all the fireflies. Moreover, this luminosity depends on the state of all the flashing phases so that it is possible to write $S_i(\phi_1, \phi_2, ..., \phi_N)$. Sensory feedback will modify the dynamics of the phase $\phi_i$ of the oscillator with the sensitivity function [32]. For an arbitrary distribution of fireflies, determining the relationship between the local luminosity $S_i$ and the set $\{I_k\}$ of emissions can be a difficult task [33].

Our experimental system features the same principle of communication. Each robot module is an elementary block of the network characterized by an oscillator of phase $\phi_i$, an action $\tau_i$ and a sensory input $S_i$. The action is provided by the torque $\tau_i$ exerted by a servomotor on the i-th link, whose phase $\phi_i$ is controlled by an oscillator with $\tau_i \propto \cos(\phi_i)$. In turn, the phase $\phi_i$ is modulated by the differential hydrodynamic forces $S_i$ resulting from the dynamics of the activation of all the servomotors $\{\tau_k\}$. In fact, an action produced by the k-th servomotor could *a priori* contribute to the motion of the i-th module, modify the local hydrodynamic load, and thus modulate the phase of the i-th oscillator.

The key difference in this analogy is the medium of communication. In fireflies, the information on the state of the oscillators is determined by the net luminosity resulting from the coherence of flashes. Since the speed of light is larger than any other speed, it is reasonable to assume that this transfer is nearly instantaneous. In contrast, in the robotic experiment, communication between oscillators is based on the transport of momentum. Each servomotor periodically injects angular momentum into the robot structure, which is then communicated to the fluid. This information is then retrieved by each module by measuring the differential force due to the exchange of momentum with the fluid. Unlike luminosity-based communication, momentum-based information transfer can be delayed and depends on the transfer function of the fluid-structure interactions or the viscoelastic properties of the robot. In addition, the presence of body deformation modes can provide long-range interactions. A servomotor located at an antinode of a standing bending wave can efficiently
transfer momentum in the serial structure of the robot, which can then be transferred to a module close to another antinode. Consequently, the complex physics of swimming determines the nature of the couplings between the oscillators. A change in gait will therefore affect the structure of this network.

Again, establishing a simple relationship between the torques and the fluid forces is a complex task. Unfortunately, the complexity of the physical situation does not allow us to obtain an explicit relationship of these couplings. However, it is still possible to assume that the differential fluid force, like the luminosity of fireflies, is an implicit function of all the phases with $S_i(\phi_1, \phi_2, ..., \phi_N)$.

This preliminary hypothesis allows us to develop a framework for studying the stability of the network. In the next section, we show how to write the complex network of oscillators from this statement, and identify the coupling strength and delay from the generalized Fourier decomposition of $S_i(\phi_1, \phi_2, ..., \phi_N)$. In general, these couplings can lead to very complex, even chaotic, dynamics. In section VC, we show that our experimental configuration, shared by several other experiments [6, 17], simplifies these dynamics. Indeed, if the amplitude of the diffusive couplings $C$ in the CPG is several orders of magnitude larger than sensory feedback, characterized by the coefficient $\sigma$, then the network dynamics can be linearized. In this case, we observe that these leading order couplings produce only a frequency shift of the oscillators, while the phase entrainment phenomenon is negligible. This result is in agreement with our experimental data and explains the robustness of this scheme.

## B. Mathematical formulation of the complex network

In this section, we propose to establish mathematically the origin of the complex network of oscillators produced by the feedback by using the slow phase approximation [22, 34] to identify the dominant terms in the phase dynamics. We assume that the hydrodynamic load $S_i(t)$ acting on the i-th module is an implicit function of the state of the phases given by $(\phi_1, \phi_2, ..., \phi_N)$. Although it is difficult to express this function explicitly, we can exploit the periodicity of the phases so that $S_i$ can be represented by a $2\pi$-periodic function of all the phases $\phi_k$ with $k \in [1, 10]$ reading in complex form [34]
$$S_i(\phi_1, \dots) = \sum_{k=1}^{10} c_{ik} e^{j\phi_k} + \sum_{k,l=1}^{10} d_{ikl} e^{j\phi_k + j\phi_l} + \text{c.c.} + \dots, \quad (10)$$

where $j$ is the unitary imaginary number. This equation is an extension of the usual Fourier series, in which the dynamics is decomposed as a function of the phases $\phi_i(t)$, which govern the periodicity of the physical signals. The term c.c. stands for the complex conjugate. We have written only the first higher harmonics of the Fourier series whose angular frequency is $2\Omega$, while the dots represent the remaining higher harmonics. Note that each Fourier coefficient $c_{ik} \in \mathbb{C}$ can be interpreted as the linear transfer function between the torque $\tau_k \propto \cos(\phi_k)$ on the $k$th joint and the differential force $S_i$ acting on the $i$th segment. These coefficients $c_{ik}$ will determine the structure of the oscillator network.

The modulation of the phase dynamics by the feedback in Eq. (2) is given by the product of the sensitivity function $Z_i = \sigma \cos(\phi_i)$ with the local feedback $S_i$, which is in complex form

$$\sigma \cos(\phi_i) S_i = \frac{\sigma}{2} \sum_{k=1}^{10} c_{ik} e^{j(\phi_k - \phi_i)} + \frac{\sigma}{2} \sum_{k=1}^{10} c_{ik} e^{j(\phi_k + \phi_i)} \\ + \frac{\sigma}{2} \sum_{k,l=1}^{10} d_{ikl} e^{j(\phi_k + \phi_l - \phi_i)} + \dots \\ + \text{c.c.} \qquad (11)$$

When $\sigma$ is small but not zero, we expect a deviation of the phase difference $(\phi_k - \phi_i)$ from zero, i.e. from the fully synchronized state with $\phi_i = \omega_0 t$, for all $i \in [1, N]$. This deviation occurs on a slow time scale $1/\sigma$ [34] with $\Omega - \omega_0 = O(\sigma)$, where $\Omega$ is the mutual frequency and $\omega_0$ the intrinsic frequency. In contrast, the second and third terms of Eq. (11) vary in comparison on fast time scales given by $2\pi/2\omega_0$ and $2\pi/\omega_0$, respectively. The standard phase reduction [22, 34] shows that the contribution of fast terms leads to a phase deviation of order of $O(\sigma)$ with respect to the slow terms. To preserve the essential dynamics [34], the fast terms are removed by averaging equation (11), which is consistent with our data analysis.

Finally, we obtain the complex network of oscillators by averaging the correlation $Z_i S_i$ given by Eq. (11) and introducing it into Eq. (2). The resulting slow phase dynamics can be then written as
$$ \dot{\phi}_i = \omega_0 + C \sum_{k=\pm 1} \sin(\phi_{i+k} - \phi_i) \\ + \sigma \sum_{k=1}^{N} |c_{ik}| \cos[\phi_k - \phi_i + \arg(c_{ik})] + O(\sigma^2). \qquad (12) $$

We obtain a network with symmetric couplings of amplitude C, which prevent the appearance of a traveling wave in the network, superimposed on nonuniform couplings with weights $\sigma|c_{ik}|$ and delays $\arg(c_{ik})$, which can be responsible for the transition. This type of complex network of oscillators can describe instabilities [10], and we assume that a change in the network couplings is at the origin of the sudden gait transition. Although the stability analysis can be performed by various methods [10, 24], they require explicit knowledge of the complex coefficients $c_{ik} \in \mathbb{C}^{N \times N}$. Unfortunately, we only have access to the N phase shifts $\{\phi_i\}$ to compute the $2N^2$ real coupling coefficients that define the real and imaginary components of $c_{ik}$. Instead of solving this complex inverse problem, we will take advantage of the smallness of the parameter $\sigma$ with respect to C to analyze the network close to the threshold $\sigma_c$.

### C. Pertubative approach for the network

Our approach is to expand the phase shift $\delta\psi_i$ and the frequency shift $\delta\Omega = \Omega - \omega_0$ of the phase $\phi_i$, as an asymptotic series with small parameter $\epsilon = \sigma/C$. This parameter quantifies the relative strength of the network perturbation by the feedback with respect to the stabilizing diffusive scheme. Note that this parameter remains small in our setup with $\epsilon = O(10^{-2})$. Thus, the mutual frequency $\Omega$ and the phase lag $\delta\psi_i$ are expanded as a power series of $\epsilon$ given by

$$ \Omega = \omega_0 + \epsilon\delta\Omega^{(1)} + \epsilon^2\delta\Omega^{(2)} + O(\epsilon^2), \\ \delta\psi_i(t) = \epsilon\delta\psi_i^{(1)} + \epsilon^2\delta\psi_i^{(2)} + O(\epsilon^2), \qquad (13) $$

with $\delta\Omega^{(k)}$ and $\delta\psi_i^{(k)}$ the kth-order correction. The synchronized state without phase lag is recovered for $\epsilon = 0$. We also keep the time dependence of the phase shift to study the stability of the network. Introducing this perturbative development in the complex network, we obtain the following leading order approximation for the averaged correlation $\overline{Z_i S_i}$
$$ \overline{Z_i S_i} = \delta\omega_i + \epsilon \sum_{k=1}^{N} L_{ik} (\delta\psi_k^{(1)} - \delta\psi_i^{(1)}) + O(\epsilon), \quad (14) $$

where the dimensionless frequency detuning $\delta\omega_i$ reads

$$ \delta\omega_i = \sum_{k=1}^{N} |c_{ik}| \cos(\arg(c_{ik})). \quad (15) $$

Similarly, the coupling coefficients are given by

$$ L_{ik} = -|c_{ik}| \sin(\arg(c_{ik})). \quad (16) $$

The results show that the leading effect of sensory feedback is to detune the intrinsic eigenfrequency by $\delta\omega_i$ for each oscillator. In contrast, the linear coupling terms producing phase entrainment are only relevant at order $\epsilon^2$, which is of the same order as the fast terms identified in the previous section, and can be neglected at leading order.

Before introducing Eq. (14) into averaged phase equation (12), the problem is nondimensionalized via the change of variable $t \to \tau = Ct$, as suggested by the linear analysis for $\sigma = 0$ (cf section III B). At leading order, we obtain a linear equation for the phase shift $\delta\psi_i^{(1)}$ and the frequency shift $\delta\Omega^{(1)}$ by introducing the phase expansion (13) into the slow phase dynamics (cf Eq.(12)) and the averaged correlation to obtain the following system

$$ \frac{d}{d\tau}\delta\psi_i^{(1)} + \delta\Omega^{(1)} = (\delta\psi_{i+1}^{(1)} + \delta\psi_{i-1}^{(1)} - 2\delta\psi_i^{(1)}) + \delta\omega_i. \quad (17) $$

This equation remains linearly stable with respect to $\delta\psi_i^{(1)}$ since its linear part corresponds to the diffusive coupling already studied in section III B. Using the eigenvalues $\hat{C}_n$ and eigenmodes $\mathbf{v}_n$ of the diffusive scheme, the final phase lag $\delta\psi_i^{(1)}$ and mutual frequency increase $\delta\Omega_i^{(1)}$ can be calculated with

$$ \delta\Omega^{(1)} = \frac{\epsilon}{N} \sum_{i=1}^{N} \delta\omega_i, \quad \delta\hat{\psi}_n^{(1)} = \frac{\epsilon}{|\hat{C}_n|} \langle \mathbf{v}_n, \delta\omega \rangle. \quad (18) $$

Here, $\delta\omega$ is the vector with component $\delta\omega_i$. Finally, the local phase $\phi_i$ after synchronization can be written into a dimensionalize form as

$$ \phi_i = (\omega_0 + \sigma\delta\Omega^{(1)}) t + \sigma \sum_n \delta\hat{\psi}_n^{(1)} v_n(i) + O(\sigma). \quad (19) $$
Consequently, the solution of the phase equation for moderate $\sigma$ can be computed from the frequency shift $\delta\omega_i$ along the oscillator chain. Note that the estimation of $\delta\omega_i$ doesn't require the explicit knowledge of the values of the coupling coefficients $c_{ik}$. Instead of using Eq. (15), the frequency shift can be obtained directly via Eq. (14) using the dominant term of the averaged correlation with $\overline{Z_i S_i} \simeq \delta\omega_i$. In general, this term is neglected in order to directly evaluate the effect of the phase entrainment associated with the $L_{ik}$ coefficients on the stability of the network [23, 34]. However, we will show that this term is sufficient to explain the observed frequencies and phases.

We have now shown that the final state of the oscillator network, i.e. the mutual frequency and phase shift, can be estimated after synchronization using the intrinsic frequency detuning $\delta\omega_i$. Compared to previous studies, this simplified model provides a better understanding of how sensory feedback works by explicitly linking the perturbation of the CPG given by $Z_i S_i$ to the state of the system. In the limit $\sigma/C$ small, the correlation $\overline{Z_i S_i}$ leads to a detuning $\delta\omega_i$ of the intrinsic frequency, which acts as a source term in a heat equation on a discrete domain defined by the diffusive couplings of the CPG.

### D. Comparison with the experimental data

To compute the intrinsic frequency detuning $\delta\omega_i$ from the experimental data, we start with the result given in equation 14. This result shows that the correlation $\overline{Z_i S_i}$ is at leading order equal to $\delta\omega_i$. Therefore, for each sensor, we calculate the time average of $Z_i S_i$ and assign its value to $\delta\omega_i$. The validity of this truncation will be justified a posteriori.

We have reported in Fig. 13 the intrinsic frequency detuning $\sigma\delta\omega_i$ along the chain of oscillators for different $\sigma$. We note in Fig. 13 that the frequency shift $\delta\omega_i$ is of order $O(1)$ for $\sigma = 0.1$, i.e. before instability, which confirms that the correlation $\overline{Z_i S_i}$ does not vanish even for small phase shifts $\psi_{i+k} - \psi_i$. One also observes the formation of a gradient of the frequency shift along the oscillator chain as $\sigma$ increases. Besides, for a given oscillator, $\sigma\delta\omega_i$ is a positive and monotonically increasing function of $\sigma$, except for $i = 10$.

Based on the results of the previous section, one can compute the modal amplitudes $\delta\hat{\psi}_n$ and net the frequency variation $\delta\Omega$ from the frequency shift $\delta\omega_i$ via Eq. (18). We have reported in Fig. 6 the estimated frequency variation $\delta\Omega$ (black crosses), and the modal amplitudes $\delta\hat{\psi}_1$ and $\delta\hat{\psi}_2$ in Fig. 7 (b) (blue and red crosses, respectively). We observe a
<figure><img src="image_15.png" /><figcaption>FIG. 13. Frequency detuning $\sigma\delta\omega_i$ for different $\sigma$ computed from $\overline{Z_iS_i}$. These curves are associated with the curves $\delta\psi_i$ in Fig. 7.(a). The dotted curve corresponds to $\sigma = 0.1$ below the onset, while the thick black curves correspond to $\sigma = [0.18, 0.21, 0.3, 0.5, 1]$.</figcaption></figure>

very good agreement between the estimated and the measured values. Surprisingly, this agreement remains valid far beyond the threshold of the instability $\sigma_c$. These results show that our analysis allows us to calculate the resulting phase shifts for $\sigma/C$ small by knowing only the average correlations $\overline{Z_iS_i}$ without explicit knowledge of the coupling coefficients.

The increase of the frequency variation $\delta\Omega$ in Fig. 6 results from the spatial average of the frequency shift $\delta\omega_i$, and is thus a direct consequence of the monotonic increase of each frequency shift $\delta\omega_i$. Similarly, the phase gradient $\delta\psi_i$ at the origin of the undulatory swimming (see Fig. 7.(a)) results from the gradient of the frequency shift, whose spatial small-scale fluctuations are smoothed by the diffusive couplings.

This result confirms that sensory feedback leads to a frequency detuning of the intrinsic frequencies of the oscillators and acts as a source term in the discrete diffusion equation for the phase. Indeed, if the phase entrainment characterized by the couplings $L_{ik}$ in Eq. (14) were of the same order of magnitude as those of the diffusive scheme, then we would not have been able to estimate $\delta\omega_i$, and the components $\delta\hat{\psi}_n$ and $\delta\Omega$. Moreover, this result also rules out the possibility of sensory feedback $S_i$ independent of the phase dynamics, a scenario previously proposed [30], as it would have produced a correlation $Z_iS_i$ dependent on the local phase dynamics $\phi_i$. These results prove that there is no phase entrainment due to sensory feedback.

This paper combines experimental and theoretical evidence for frequency self-tuning in a CPG network via sensory feedback at the origin of gait selection.
## E. Stability of the body phase dynamics

We have demonstrated that the observed instability in the network cannot be explained by a mechanism of phase entrainment, since the phase equation remains stable due to diffusive coupling. If the oscillator network and the body dynamics are stable when both systems are isolated, the origin of the instability is very likely due to the couplings between these two subsystems as hypothesized in the introduction. In this section, we propose to explore this hypothesis by deriving a model that accounts for the instability and is consistent with our experimental data.

First, we observed that the growth rate $\xi$ of the instability (see Fig. 9) does not vary linearly with $\sigma - \sigma_c$ but describes a nonlinear behavior with a square root-like dependence. This property excludes standard one-dimensional supercritical bifurcations (such as the pitchfork bifurcation). Indeed, the Taylor expansion of the growth rate $\xi$ of a supercritical bifurcation would result in $\xi \propto \sigma - \sigma_c$, which is not in agreement with the observations.

Among all the possible instability scenarios, we propose to extend the state variable to two dimensions by considering the odd modes that break the symmetry of the robot dynamics: one for the phase lag $\delta\hat{\psi}_1$ and the other one for the body deformation $|\hat{\theta}_1|$. The linear evolution of the state variable $\mathbf{x}(t) \in \mathbb{R}^2$ with $\mathbf{x}(t) = (|\hat{\theta}_1|, \delta\hat{\psi}_1)$ is thus given by

$$ \dot{\mathbf{x}}(t) = A(\sigma)\mathbf{x}, \qquad (20) $$

with $A \in \mathbb{R}^{2 \times 2}$ a matrix with real coefficients $A_{ik}$. Due to the complexity of the nonlinear experimental system, we cannot explicitly identify the values and the physical mechanisms at the origin of these couplings. A new study based on numerical simulations could facilitate the identification of these couplings by having access to all the physical quantities of the problem.

Nevertheless, it is still possible to study the implications of such a minimal linear model that couples phase and body dynamics. First, we assume that each mode is linearly stable in the absence of coupling, so that the terms $A_{ii}$ are negative. For convenience, we rewrite these terms as $-|A_{ii}|$. Second, we state that motor activation controlled by the oscillator network leads to the growth of the odd mode $|\hat{\theta}_1|$. We assume here that it takes the form of a forcing term given by $A_{12}\delta\hat{\psi}_1$, so that an odd component in the motor coordination leads to an odd component in the body deformation. Finally, we hypothesize that a symmetry
breaking in the body dynamics will be reflected in the distribution of hydrodynamic forces along the body. This symmetry breaking should also be present in the frequency shift $\delta\omega_i$ produced by the feedback. This hypothesis leads to a linear relation between $\delta\omega_i$ and $|\hat{\theta}_1|$, which takes the form: $\langle \mathbf{v}_1, \delta\omega \rangle = A_{21}|\hat{\theta}_1|$ where $A_{21}$ is also constant.

Finally, system 20 and the matrix A can be written

$$ \frac{d}{dt} \begin{pmatrix} |\hat{\theta}_1| \\ \delta\hat{\psi}_1 \end{pmatrix} = \begin{pmatrix} -|A_{11}| & A_{12} \\ \sigma A_{21} & -|A_{22}| \end{pmatrix} \begin{pmatrix} |\hat{\theta}_1| \\ \delta\hat{\psi}_1 \end{pmatrix}. \quad (21) $$

Note the presence of the prefactor $\sigma$ that accounts for the strength of sensory feedback in the phase dynamics. For $\sigma = 0$, the system is stable such that the mode $(\hat{\theta}_1, \hat{\psi}_1)$ converges to zero, as is the case in our experiment if we neglect the small symmetry imperfection. This linear system can become unstable for increasing $\sigma$ while keeping the coefficients $A_{ik}$ constant, if one of the eigenvalues of the matrix becomes positive. This case occurs when $A_{12}A_{21}$ is positive leading to positive feedback. To overcome the damping of the modes, the feedback strength must be $\sigma > \sigma_c$ with

$$ \sigma_c = \frac{|A_{11}| |A_{22}|}{A_{12} A_{21}}. \quad (22) $$

In this condition, the growth rate of the instability can be written in a simple form given by

$$ \xi = -\nu \pm \sqrt{\nu^2 + \alpha(\sigma - \sigma_c)}, \quad (23) $$

with the coefficients $\nu$ and $\alpha$ defined by

$$ \nu = \frac{|A_{11}| + |A_{22}|}{2}, \quad \alpha = A_{12}A_{21}. \quad (24) $$

We show in Fig. 9 that the grow rate given by Eq. (23) reproduces well the behavior of the observed growth rate by fixing the free parameters to $\nu = 0.0165 \text{ s}^{-1}$ and $\alpha = 0.9 \text{ s}^{-1}$.

This simple toy model allows us to explain how two initially stable systems given by the body and the CPG can become unstable by increasing the strength of sensory feedback, as in the experiment. This model also reproduces the behavior of the observed growth rate $\xi$, which does not vary linearly with the parameter $\sigma$ close to the threshold $\sigma_c$. Note that the hypothesis of constant linear coupling coefficients is only valid close to the onset of the
instability and for small amplitudes of $\delta\hat{\psi}_1$ and $\hat{\theta}_1$. In general, these coupling coefficients $A_{ik}$ are strongly dependent on gait properties and therefore vary with $\sigma$.

This model does not exclude other scenarios, and requires further studies to be validated. However, it illustrates well how a symmetry breaking in the neural system can induce a symmetry breaking in the body dynamics and vice versa, as shown in Fig. 2.

# VI. CONCLUSION

In this study, we identified a number of mechanisms explaining the impact of sensory feedback on the central pattern generator, such as oscillator frequency detuning, large-scale phase dynamics and reciprocal coupling of modes in the oscillators and body dynamics.

Until now, it was difficult to detect these elements because they were masked by the non-symmetrical CPG, which naturally produces phase shifts along the spinal cord. The implementation of a symmetrical CPG has therefore proved fruitful and could be extended to other types of animal locomotion or other environments (granular or viscous) via robotic or numerical platforms.

Our theoretical study motivated by the experimental results has also shown that sensory feedback produces at leading order a steady frequency detuning of the oscillator along the spinal cord. Surprisingly, such mechanism can also be naturally present in other biological systems [35, 36]. In lampreys, proprioceptive feedback can cause a *tonic excitatory effect*, leading to an increase of the mutual frequency in the network [37, 38]. In our study, we have therefore identified a possible source of this tonic excitatory effect which in our experiment comes from an average shift of all the intrinsic frequencies (see section V). Our analytical study can also justify the *adhoc* introduction of a tonic excitatory effect in previous numerical studies [7, 14].

We believe that the mechanisms identified will still be present for more realistic configurations of the CPGs, even if other phenomena may be superimposed on them. A natural continuation of this study would be to analyze these solutions by progressively implementing an intrinsic phase gradient in the CPGs, for example, thanks to a phase shift [17]. Moreover, it will be interesting to derive explicitly the toy model and better characterize the interaction between the body and phase dynamics.
## VII. ACKNOWLEDGMENT

The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 730994 (TER-RINet). This work was also supported by Human Frontier Science Program (HFSP), grant RGP0027/2017, by the ERC synergies Grant Salamandra (945517), and by the French National Research Agency ANR (grant no. ANR-2019-CE33-0004-01).

## VIII. BILBIO

[1] M. S. Triantafyllou and G. S. Triantafyllou, An efficient swimming machine, Sci. Am. **272**, 64 (1995).

[2] A. J. Ijspeert, A. Crespi, D. Ryczko, and J.-M. Cabelguen, From swimming to walking with a salamander robot driven by a spinal cord model, Science **315**, 1416 (2007).

[3] J. Aguilar, T. Zhang, F. Qian, M. Kingsbury, B. McInroe, N. Mazouchova, C. Li, R. Maladen, C. Gong, M. Travers, et al., A review on locomotion robophysics: the study of movement at the intersection of robotics, soft matter and dynamical systems, Rep. Prog. Phys. **79**, 110001 (2016).

[4] A. Anastasiadis, L. Paez, K. Melo, E. Tytell, A. J. Ijspeert, and K. Mulleners, Identification of the trade-off between speed and efficiency in undulatory swimming using a bio-inspired robot, Scientific Reports **13**, 15032 (2023).

[5] M. Gazzola, M. Argentina, and L. Mahadevan, Gait and speed selection in slender inertial swimmers, Proc. Natl. Acad. Sci. **112**, 3874 (2015).

[6] E. D. Tytell, M. C. Leftwich, C.-Y. Hsu, B. E. Griffith, A. H. Cohen, A. J. Smits, C. Hamlet, and L. J. Fauci, Role of body stiffness in undulatory swimming: insights from robotic and computational models, Phys. Rev. Fluids **1**, 073202 (2016).

[7] C. Hamlet, L. Fauci, J. R. Morgan, and E. D. Tytell, Proprioceptive feedback amplification restores effective locomotion in a neuromechanical model of lampreys with spinal injuries, Proceedings of the National Academy of Sciences **120**, e2213302120 (2023).
[8] J. Collins and I. Stewart, Hexapodal gaits and coupled nonlinear oscillator models, Biological cybernetics **68**, 287 (1993).

[9] H. Abarbanel, M. I. Rabinovich, A. Selverston, M. Bazhenov, R. Huerta, M. Sushchik, and L. Rubchinskii, Synchronisation in neural networks, Physics-Uspekhi **39**, 337 (1996).

[10] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Synchronization in complex networks, Phys. Rep. **469**, 93 (2008).

[11] A. J. Ijspeert, Central pattern generators for locomotion control in animals and robots: a review, Neur. Net. **21**, 642 (2008).

[12] S. Grillner, The motor infrastructure: from ion channels to neuronal networks, Nature Reviews Neuroscience **4**, 573 (2003).

[13] A. H. Cohen and P. Wallén, The neuronal correlate of locomotion in fish, Exp. Brain Res. **41**, 11 (1980).

[14] C. L. Hamlet, K. A. Hoffman, E. D. Tytell, and L. J. Fauci, The role of curvature feedback in the energetics and dynamics of lamprey swimming: A closed-loop model, PLoS computational biology **14**, e1006324 (2018).

[15] A. Frigon and S. Rossignol, Experiments and models of sensorimotor interactions during locomotion, Biological cybernetics **95**, 607 (2006).

[16] J. Knüsel, A. Crespi, J.-M. Cabelguen, A. J. Ijspeert, and D. Ryczko, Reproducing five motor behaviors in a salamander robot with virtual muscles and a distributed cpg controller regulated by drive signals and proprioceptive feedback, Frontiers in Neurorobotics **14**, 604426 (2020).

[17] R. Thandiackal, K. Melo, L. Paez, J. Herault, T. Kano, K. Akiyama, F. Boyer, D. Ryczko, A. Ishiguro, and A. J. Ijspeert, Emergence of robust self-organized undulatory swimming based on local hydrodynamic force sensing, Science Rob. **6** (2021).

[18] J. Sánchez-Rodríguez, F. Celestini, C. Raufaste, and M. Argentina, Proprioceptive mechanism for bioinspired fish swimming, Physical Review Letters **126**, 234501 (2021).

[19] M. C. Leftwich, E. D. Tytell, A. H. Cohen, and A. J. Smits, Wake structures behind a swimming robotic lamprey with a passively flexible tail, Journal of Experimental Biology **215**, 416 (2012).

[20] J. Christenson, A. Boman, P.-Å. Lagerbäck, and S. Grillner, The dorsal cell, one class of primary sensory neuron in the lamprey spinal cord. i. touch, pressure but no nociception—a physiological study, Brain research **440**, 1 (1988).
[21] A. H. Cohen, G. B. Ermentrout, T. Kiemel, N. Kopell, K. A. Sigvardt, and T. L. Williams, Modelling of intersegmental coordination in the lamprey central pattern generator for locomotion, *Trends in neuro. 15*, 434 (1992).

[22] Y. Kuramoto, *Chemical oscillations, waves, and turbulence* (2003).

[23] L. M. Pecora and T. L. Carroll, Master stability functions for synchronized coupled systems, *Phys. Rev. Lett.* **80**, 2109 (1998).

[24] K. S. Fink, G. Johnson, T. Carroll, D. Mar, and L. Pecora, Three coupled oscillators as a universal probe of synchronization stability in coupled oscillator arrays, *Physical Review E* **61**, 5080 (2000).

[25] A. J. Smits, Undulatory and oscillatory swimming, *J. Fluid Mech.* **874** (2019).

[26] P. Holmes, J. L. Lumley, and G. Berkooz, *Turbulence, coherent structures, symmetry and dynamical systems* (2010).

[27] G. M. L. Gladwell, The inverse problem for the vibrating beam, *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences* **393**, 277 (1984).

[28] S. Ramananarivo, R. Godoy-Diana, and B. Thiria, Propagating waves in bounded elastic media: Transition from standing waves to anguilliform kinematics, *EPL* **105**, 54003 (2014).

[29] F. Pétrélis, C. Laroche, B. Gallet, and S. Fauve, Drifting patterns as field reversals, *EPL* **112**, 54007 (2015).

[30] D. Owaki, T. Kano, K. Nagasawa, A. Tero, and A. Ishiguro, Simple robot suggests physical interlimb communication is essential for quadruped walking, *Journal of The Royal Society Interface* **10**, 20120669 (2013).

[31] A. T. Winfree, *The geometry of biological time*, Vol. 2 (Springer, 1980).

[32] S. H. Strogatz, From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators, *Physica D* **143**, 1 (2000).

[33] M. McCrea, B. Ermentrout, and J. E. Rubin, A model for the collective synchronization of flashing in photinus carolinus, *Journal of the Royal Society Interface* **19**, 20220439 (2022).

[34] A. Pikovsky, M. Rosenblum, and J. Kurths, *Synchronization: a universal concept in nonlinear science* (2002).

[35] N. E. Diamant and A. Bortoff, Nature of the intestinal low-wave frequency gradient, *Amer. J. Physiol.* **216**, 301 (1969).
[36] G. B. Ermentrout and N. Kopell, Frequency plateaus in a chain of weakly coupled oscillators, i., SIAM **15**, 215 (1984).

[37] L. Guan, T. Kiemel, and A. H. Cohen, Impact of movement and movement-related feedback on the lamprey central pattern generator for locomotion, Journal of Experimental Biology **204**, 2361 (2001).

[38] T. Kiemel and A. H. Cohen, Bending the lamprey spinal cord causes a slowly-decaying increase in the frequency of fictive swimming, Brain research **900**, 57 (2001).