# Reproducing Five Motor Behaviors in a Salamander Robot With Virtual Muscles and a Distributed CPG Controller Regulated by Drive Signals and Proprioceptive Feedback

Jeremie Knusel, Alessandro Crespi, Jean-Marie Cabelguen, Auke J. Ijspeert, Dimitri Ryczko

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

To cite this version:

Jeremie Knusel, Alessandro Crespi, Jean-Marie Cabelguen, Auke J. Ijspeert, Dimitri 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, 2020, 14, 10.3389/fnbot.2020.604426. [https://doi.org/10.3389/fnbot.2020.604426](https://doi.org/10.3389/fnbot.2020.604426)

HAL Id: hal-04125649

https://hal.science/hal-04125649

Submitted on 12 Jun 2023

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
<figure><img src="image_3.png" /></figure>

# Reproducing Five Motor Behaviors in a Salamander Robot With Virtual Muscles and a Distributed CPG Controller Regulated by Drive Signals and Proprioceptive Feedback

Jérémie Knüsel<fnref n="1" />,², Alessandro Crespi<fnref n="1" />, Jean-Marie Cabelguen<fnref n="3" />†, Auke J. Ijspeert<fnref n="1" />* and Dimitri Ryczko⁴,⁵,⁶,⁷*

## OPEN ACCESS

**Edited by:**

Poramate Manoonpong,
University of Southern
Denmark, Denmark

**Reviewed by:**

Cecilia Laschi,
Sant'Anna School of Advanced
Studies, Italy
Shinya Aoi,
Kyoto University, Japan
Holk Cruse,
Bielefeld University, Germany

**\*Correspondence:**

Dimitri Ryczko
dimitri.ryczko@usherbrooke.ca
Auke J. Ijspeert
auke.ijspeert@epfl.ch

Received: 09 September 2020

Accepted: 10 November 2020

Published: 23 December 2020

**Citation:**

Knüsel J, Crespi A, Cabelguen J-M, Ijspeert AJ and Ryczko D (2020)

Reproducing Five Motor Behaviors in a Salamander Robot With Virtual Muscles and a Distributed CPG Controller Regulated by Drive Signals and Proprioceptive Feedback.

Front. Neurorobot. 14:604426.

doi: 10.3389/fnbot.2020.604426

Diverse locomotor behaviors emerge from the interactions between the spinal central pattern generator (CPG), descending brain signals and sensory feedback. Salamander motor behaviors include swimming, struggling, forward underwater stepping, and forward and backward terrestrial stepping. Electromyographic and kinematic recordings of the trunk show that each of these five behaviors is characterized by specific patterns of muscle activation and body curvature. Electrophysiological recordings in isolated spinal cords show even more diverse patterns of activity. Using numerical modeling and robotics, we explored the mechanisms through which descending brain signals and proprioceptive feedback could take advantage of the flexibility of the spinal CPG to generate different motor patterns. Adapting a previous CPG model based on abstract oscillators, we propose a model that reproduces the features of spinal cord recordings: the diversity of motor patterns, the correlation between phase lags and cycle frequencies, and the spontaneous switches between slow and fast rhythms. The five salamander behaviors were reproduced by connecting the CPG model to a mechanical simulation of the salamander with virtual muscles and local proprioceptive feedback. The main results were validated on a robot. A distributed controller was used to obtain the fast control loops necessary for implementing the virtual muscles. The distributed control is demonstrated in an experiment where the robot splits into multiple functional parts. The five salamander behaviors were emulated by regulating the CPG with two descending drives. Reproducing the kinematics of backward stepping and struggling however required stronger muscle contractions. The passive oscillations observed in the salamander's tail during forward underwater stepping could be reproduced using a third descending drive of zero to the tail oscillators. This reduced the drag on the body in our hydrodynamic simulation. We explored the effect of local proprioceptive feedback

<footnote n="1">¹Biorobotics Laboratory (BioRob), Institute of Bioengineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland, ²Institute for Optimisation and Data Analysis (IODA), Bern University of Applied Sciences, Biel, Switzerland,</footnote>

<footnote n="3">³Institut National de la Santé et de la Recherche Médicale (INSERM) U 862 - Neurocentre Magendie, Université de Bordeaux, Bordeaux, France, ⁴Département de Pharmacologie-Physiologie, Faculté de Médecine et des Sciences de la Santé, Université de Sherbrooke, Sherbrooke, QC, Canada, ⁵Centre de Recherche du Centre Hospitalier Universitaire de Sherbrooke, Sherbrooke, QC, Canada, ⁶Institut de Pharmacologie de Sherbrooke, Sherbrooke, QC, Canada, ⁷Centre d'Excellence en Neurosciences de l'Université de Sherbrooke, Sherbrooke, QC, Canada</footnote>

<footnote>†Deceased</footnote>
during swimming and forward terrestrial stepping. We found that feedback could replace or reduce the need for different drives in both cases. It also reduced the variability of intersegmental phase lags toward values appropriate for locomotion. Our work suggests that different motor behaviors do not require different CPG circuits: a single circuit can produce various behaviors when modulated by descending drive and sensory feedback.

**Keywords:** central pattern generator (CPG), proprioceptive sensory feedback, descending drive, distributed control, salamander, locomotion, numerical modeling, robotics

# INTRODUCTION

Many motor behaviors in animals require coordinated rhythmic activation of multiple muscles. Neural networks capable of producing such activity patterns without rhythmic input from other networks or from sensory feedback are called central pattern generators (CPGs). It has been shown that CPGs in the spinal cord underlie locomotion in many vertebrate species (for review, Grillner and El Manira, 2020). Drive signals descending from brain neurons control locomotion initiation, speed and gait transitions (Brocard et al., 2010; Capelli et al., 2017; Caggiano et al., 2018; Josset et al., 2018) and steering movements (Fagerstedt et al., 2001; Ryczko et al., 2016b; Cregg et al., 2020). Sensory feedback plays an important role in modulating the CPG activity to adapt the locomotor pattern to the environment (e.g., Wyart et al., 2009; Akay et al., 2014; Hubbard et al., 2016; Knafo et al., 2017). These feedback signals depend on the interactions between the neural networks, the mechanical properties of the body and the environment, making it a challenge to fully understand the operation of the CPG even at a high level of abstraction. Numerical models of the complete system can be used to investigate the effect of sensory feedback on the CPG, but some aspects such as hydrodynamic and friction forces are difficult to simulate reliably. Robots are thus useful to validate simulation results in the real world, with real physics. Here, we used numerical simulations and robotics to investigate the generation of different behaviors in the salamander, an interesting animal model as it can move underwater and on ground (Ryczko et al., 2020). In particular, we addressed the following questions:

1. Can different motor behaviors be generated by a single spinal CPG circuit as opposed to requiring several dedicated CPG circuits?

2. What are the roles of descending drives in generating these different motor behaviors, and how many independent drives are necessary?

3. What is the potential role of sensory feedback in shaping the patterns and in reducing the variability of CPG activity observed in isolated spinal cords?

We used the *Salamandra robotica II* robot (Crespi et al., 2013) driven by a spinal CPG model and virtual muscles to reproduce the five salamander behaviors documented in Ryczko et al. (2015): forward swimming, forward and backward terrestrial stepping, forward underwater

stepping, and struggling. To match the biological data from that study, we focused on reproducing the patterns of muscle activation and body curvature along the body axis.

For the CPG, our starting point was the abstract oscillator model of Ijspeert et al. (2007), with modifications to allow for the flexible coordination of limb and axial network activities (Knüsel et al., 2013). This flexibility is required to reproduce the observations from Ryczko et al. (2015). Here, we extended the model to comprise 25 segments and introduced random parameters to account for the differences between individuals. The main hypotheses are the following: (1) limb oscillators project only to the axial oscillators close to the corresponding girdles; (2) couplings between axial oscillators are stronger in the head-to-tail direction; (3) limb oscillators saturate<fnref n="1" /> at lower excitatory drives than axial oscillators; (4) hindlimb oscillators are intrinsically slower than forelimb oscillators. Hypotheses 1 and 2 make the model’s intersegmental phase lag flexible and controllable (Knüsel et al., 2013). Hypotheses 3 and 4 allow the model to reproduce the distribution of phase lags of recordings *in vitro*.

We modeled the biomechanical properties of the body axis using virtual muscles that determine the torques of the axial joints based on the CPG activity and the current joint position and velocity (Ekeberg, 1993). The joint positions were also used for proprioceptive feedback, simulating stretch receptors that send phasic inputs to the local CPG segments. The virtual muscle model requires a small time step for stability and accuracy of the numerical integration, which is challenging to achieve with eight joints given the limited bandwidth and processing power of the robotic platform. We solved this difficulty by distributing the computation of the CPG and muscle models in the eight active modules, so that each module calculates the part of the model that controls its own joint. The modules use peer-to-peer communication, such that splitting the robot results in several functional pieces (unlike most robots).

The isolated CPG model was tuned to reproduce the diversity of coordination patterns observed in isolated spinal cords. To determine values for the proprioceptive feedback and virtual muscle parameters, we systematically explored the parameter space using a mechanical simulation of the robot: Using the same tonic drive for all oscillators, we identified parameter values

<footnote n="1">¹When the excitatory inputs to the spinal cord increase beyond a certain level, CPG units saturate: they cease to produce rhythmic output. In the model, this is implemented as a decrease of oscillation amplitude toward zero (see the left part of Figure 3 for an illustration).</footnote>
that allowed feedback to have a positive effect on swimming speed while maintaining stable CPG rhythms. We then attempted to reproduce the five salamander behaviors without sensory feedback, by varying the tonic drive sent to different parts of the CPG, adapting other model parameters when necessary. Further simulations were made to investigate the effect of proprioceptive feedback during swimming and forward terrestrial stepping. Finally, we used the real robot to validate simulation results. Each behavior was reproduced on the robot using different “individuals” from the family of models used in modeling isolated spinal cords, to check the robustness of the control architecture to individual variations.

## MOTOR CONTROL IN SALAMANDERS

The salamander spinal CPG produces the rhythmic movements of the limbs (Cheng et al., 1998; Lavrov and Cheng, 2004; Ijspeert et al., 2007), trunk (Delvolvé et al., 1999; Branchereau et al., 2000; Ryczko et al., 2010), and tail (Charrier and Cabelguen, 2013). This spinal circuitry is controlled by the brainstem in salamanders as in other vertebrates (for review, see Ryczko and Dubuc, 2013). Stimulation of the salamander mesencephalic locomotor region elicits stepping at low stimulation intensities, whereas swimming requires higher intensities (Cabelguen et al., 2003). These descending commands are carried to the spinal cord by reticulospinal neurons (Ryczko et al., 2016a, see also Ryczko et al., 2020 for a recent review).

The coordination of muscles along the body axis plays an important role in salamander locomotion, to generate thrust during swimming and to maximize the stride length during terrestrial stepping (Delvolvé et al., 1997). So far, at least five salamander motor behaviors have been characterized: forward swimming, forward and backward terrestrial stepping, forward underwater stepping, and struggling (Ryczko et al., 2015). Forward terrestrial stepping generally takes the form of a walking trot, but lateral sequence walks have also been observed (reviewed in Chevallier et al., 2008, see also Ashley-Ross et al., 2009). During forward underwater stepping, the salamander progresses at the bottom of water, with periods of suspension in water without ground contact. Struggling refers to the behavior of the salamander when it is firmly grasped at the pelvic girdle. Electromyographic (EMG) recordings of multiple segments in the salamander mid-trunk show that each of the five behaviors is characterized by a specific pattern of muscle activation, in terms of cycle frequencies and intersegmental phase lags: (1) rostrocaudal waves occur during forward swimming and, with lower cycle frequencies, during backward terrestrial stepping; (2) slow caudorostral waves occur during struggling; (3) standing waves are stable during forward terrestrial stepping but more variable during forward underwater stepping (Ryczko et al., 2015).

Kinematic recordings show similar patterns of trunk curvature. However, kinematic intersegmental phase lags are significantly larger during forward terrestrial stepping and swimming (Frolich and Biewener, 1992; Ryczko et al., 2015). In other words, the delay between muscle activation and body

bending gets larger toward the tail. This suggests that the mechanical properties of body tissues play an important role during these behaviors, as suggested by a lamprey modeling study (Tytell et al., 2010).

The increasing EMG-mechanical delay toward the tail also suggests that proprioceptive feedback might have a different effect at various points along the body axis. Salamanders are known to have sensory cells that generate proprioceptive information relative to axial movements: The skin contains mechano-sensitive Merkel cells (Scott et al., 1981; Diamond et al., 1986), and some cells in the spinal cord are morphologically similar to the mechano-sensitive “edge cells” (Schroeder and Egar, 1990) that encode body bending in lampreys (Grillner et al., 1982, 1984). They also have cerebrospinal fluid contacting neurons (Kolmer-Agduhr cells, Harper and Roberts, 1993), which are active during body bending in zebrafish (Böhm et al., 2016) and provide mechanosensory input to the swimming CPG (Wyart et al., 2009; Hubbard et al., 2016, Orts-Del'Immagine et al., 2020, see also Jalalvand et al., 2016 in lampreys). The limbs are another source of proprioceptive feedback, as they contain fibers that respond to stretch similarly to muscle spindles in other species (Bone et al., 1976).

According to *in vitro* recordings of the salamander spinal cord, the isolated CPG can generate stable patterns for the three types of axial waves (caudorostral, standing and rostrocaudal waves), with occasional switches between two wave types (Ryczko et al., 2015). The intersegmental phase lags generated by the isolated CPG cover a greater range than those observed in EMG recording (−12.6 to +12.4% of a cycle duration for recordings *in vitro*, and −4.8 to +6.4% for EMG recordings), with a distribution showing three peaks centered on −9.6, −1.0, and +6.6%. The salamander CPG thus provides a flexible ground onto which sensory feedback and descending drives could act to influence the spinal motor output.

## RELATED MODELING WORK

Previous studies have modeled the CPG components using abstract oscillators (Ijspeert et al., 2005, 2007; Knüsel et al., 2013; Yin et al., 2016), single bursting neurons (Liu et al., 2018, 2020), integrate-and-fire neurons (Ijspeert, 2001; Bem et al., 2003; Harischandra et al., 2011; Knüsel et al., 2013) and detailed networks of three compartment Hodgkin-Huxley neurons (Bicanski et al., 2013).

The mechanical body of the salamander has been modeled with varying accuracy. Many models include four joints between the girdles and a single degree of freedom (DOF) per limb (Ijspeert, 2001; Ijspeert et al., 2005, 2007; Suzuki et al., 2019a) or three DOFs per limb (Harischandra et al., 2010, 2011; Liu et al., 2018, 2020). The simplest model has one of each (Yin et al., 2016), while other models have one joint between the girdles and two DOFs per limb (Zhong et al., 2018; Suzuki et al., 2019b). Bem et al. (2003) have modeled the swimming salamander as a chain of ten links, corresponding roughly to three trunk joints and no limbs. The most accurate model has five joints between the girdles and four DOFs in each limb (Karakasiliotis et al., 2016; Horvat
and Ijspeert, 2017; Horvat et al., 2017). Mechanical properties (damping and elasticity) of the body tissues were included in the muscle models used by Ijspeert (2001; 2005) Bem et al. (2003), Harischandra et al. (2010, 2011), and Liu et al. (2018, 2020), and in the controller of Suzuki et al. (2019b).

The effect of sensory feedback on the activity of the salamander CPG has only been investigated in simulation (Bem et al., 2003; Ijspeert et al., 2005; Harischandra et al., 2011; Liu et al., 2020). The role of sensory feedback in body-limb coordination has also been investigated using controllers without CPG, both in simulations (Horvat and Ijspeert, 2017) and with a robot (Suzuki et al., 2019a).

Most studies have focused on the reproduction of forward terrestrial stepping (with a walking and/or trotting gait), swimming, transitions between these behaviors, and turning. The exceptions are the works of Karakasiliotis et al. (2016) which reproduced underwater stepping in addition to swimming and forward terrestrial stepping (though using predefined joint trajectories rather than a CPG) and Liu et al. (2018) which reproduced backward terrestrial stepping in addition to forward terrestrial stepping (using dedicated networks for each gait).

**Table 1** summarizes the particularities of past studies and how they compare to the present one. To our knowledge, the present

<table><caption>TABLE 1 | Related studies.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">CPG</th><th rowspan="2">Trunk joints</th><th rowspan="2">Limb DOFs</th><th rowspan="2">Robot</th><th rowspan="2">Biomech</th><th colspan="2">Proprio. feedback</th><th colspan="2">Force feedback</th><th colspan="3">Behaviors</th><th rowspan="2">Turning</th><th rowspan="2">Transition</th></tr><tr><th>Trunk</th><th>Limbs</th><th>Trunk</th><th>Limbs</th><th>Swim</th><th>Walk</th><th>Trot</th><th>Back step</th><th>U.w. step</th><th>Struggle</th></tr></thead><tbody><tr><td>Bem et al. (2003)</td><td>IF</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Bicanski et al. (2013)</td><td>HH</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Knüsel et al. (2013)</td><td>AO+IF</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ijspeert (2001)</td><td>IF</td><td>●●●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ijspeert et al. (2005)</td><td>AO</td><td>●●●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ijspeert et al. (2007)</td><td>AO</td><td>●●●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Harischandra et al. (2010)</td><td></td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Harischandra et al. (2011)</td><td>IF</td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Yin et al. (2016)</td><td>AO</td><td>●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Karakasiliotis et al. (2016)</td><td></td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Horvat et al. (2017)</td><td></td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Horvat and Ijspeert (2017)</td><td></td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Liu et al. (2018)</td><td>BN</td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhong et al. (2018)</td><td></td><td>●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Suzuki et al. (2019a)</td><td></td><td>●●●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Suzuki et al. (2019b)</td><td>AO</td><td>●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Liu et al. (2020)</td><td>BN</td><td>●●●</td><td>●●●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>This study</td><td>AO</td><td>●●●</td><td>●</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

CPG: IF, Integrate and fire; HH, Hodgkin-Huxley; AO, Abstract oscillators; BN, Bursting neurons. Trunk joints: counting joints between girdles (a joint on the girdle counts as a half). DOFs: degrees of freedom. The number of black dots represents the number of trunk joints and limb degrees of freedom respectively. Biomech: mechanical properties from body tissues such as muscles. Walk: forward terrestrial stepping with lateral sequence walking. Trot: forward terrestrial stepping with walking trot. Colors indicate different groups of model features: mechanical model (green), sensory feedback modalities (orange), behaviors (gray).
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** | Robot with distributed controller with the spinal central pattern generator (CPG) model, axial proprioceptive feedback, descending drives and virtual muscles. (A) The robot *Salamandra robotica II.* (B) The axial CPG was divided in 8 groups (gray rectangles) to distribute the computations in the 8 robot modules with active joints. Left: tonic descending drives are applied to limb (blue), trunk (green) and tail (yellow) oscillators. Virtual stretch receptors (orange triangles) project to the 3 nearest segments with opposite ipsilateral (excitatory) and contralateral (inhibitory) weights. Feedback from the neck joint (dashed orange) was disabled for robot experiments (see Results). Black lines indicate bidirectional couplings between oscillators (see **Figure 2A**). Middle: 2 outputs $x_i, x_{i+25}$ of each group (purple horizontal arrows) govern left ($l$) and right ($r$) muscle activities $M_i$ from which the muscle model calculates an output torque $T_i$. Right: the torque $T_i$ is applied at each axial joint (orange circles). The joint position $\phi_i$ and velocity $\dot{\phi}_i$ are fed back (orange arrows) to the muscles. Virtual stretch receptors only receive $\phi_i$. The phases $\theta_i$ of limb oscillators (red horizontal arrows) determine the limb positions.</figcaption></figure>
work is the first to incorporate biomechanical properties and proprioceptive sensory feedback in a real salamander robot.

# MATERIALS AND METHODS

## CPG Model

The model was developed using the Codyn framework and exported to C code to run on the robot microcontrollers. Only the 25 most rostral axial segments (each comprising 1 left and 1 right hemisegments) are modeled out of the 40 segments that salamanders typically have (for review see, Chevallier et al., 2008). Four additional oscillators control the limbs. The 25 axial segments control the active part of the robot. The caudal half of the robot tail is a passive, flexible caudal fin (Figure 1A). Each axial hemisegmental oscillator and each limb oscillator is modeled as a phase oscillator with controllable amplitude, and the connections between oscillators are functions of the phase difference between sender and receiver:

$$ 
\begin{aligned}
\dot{\theta}_i &= 2\pi v_i + \sum_j r_j w_{ij} \sin(\theta_j - \theta_i - \varphi_{ij}) - \frac{s_i}{r_i} \sin \theta_i \\
\dot{r}_i &= a(R_i - r_i) + s_i \cos \theta_i \\
x_i &= r_i (1 + \cos \theta_i) \\
v_i &= d_i e_i \\
R_i &= d_i P(d_i, d_i^{th})
\end{aligned}
$$

A positive output $x_i$ (which determines the muscle activation) is calculated from the instantaneous phase $\theta_i$ and amplitude $r_i$. The intrinsic frequency $v_i$ is proportional to the oscillator excitability $e_i$ and to a drive $d_i$ that represents the excitation from descending drives. The intrinsic amplitude $R_i$ increases with increasing drive until it approaches a saturation threshold $d_i^{th}$ after which it decreases progressively to zero due to the sigmoid function $P(d, d_i^{th}) = \frac{1}{1+e^{b(d-d_i^{th})}}$ with b the saturation rate. The excitability $e_i$ determines the intrinsic frequency of a particular oscillator as a function of the external drive. The excitability of each oscillator is drawn from a Gaussian distribution with different means for forelimb, hindlimb and axial oscillators. The saturation thresholds of the forelimbs, hindlimbs and axial network are also drawn from a Gaussian distribution with different means for the axial and limb networks. The coupling from oscillator j to oscillator i is characterized by a strength $w_{ij}$ and phase bias $\varphi_{ij}$. The gain a determines the speed of convergence for the amplitude. The symbol $s_i$ represents the feedback signal from simulated stretch receptors (see below). The terms $-\frac{s_i}{r_i} \sin \theta_i$ and $s_i \cos \theta_i$ are the polar coordinate equivalent of adding $s_i$ to the derivative $\dot{x}$ of an oscillator in Cartesian coordinates (see Supplementary Materials for the derivation).

The network connectivity is described in Figure 2A and Table 2. Other parameter values are provided in Tables 3, 4.

## Simulations of the Isolated CPG

For simulations of *in vitro* electrophysiological recordings of the isolated salamander spinal cord reported in Ryczko et al. (2015), the same drive $d_i = d$ was used for all oscillators

to represent a tonic pharmacological stimulation, with small fluctuations over time added in the form of a mean reverting random walk: $\dot{d} = c(d_0 - d) \pm \sigma$ with $d_0$ the drive picked from a Gaussian distribution, c a convergence factor, and $\pm \sigma$ a random process yielding positive and negative steps with equal probability. Multiple simulations were performed with different random seeds to reflect the diversity of coordination patterns observed in individual spinal cord preparations.

## Muscle Model

A linear spring-damper model with variable stiffness (Ekeberg, 1993) was used to model a pair of antagonist muscles and calculate the resulting torque at each axial joint (Figure 1B):

$$ T_i = \alpha (M_i^l - M_i^r) - \beta (M_i^l + M_i^r + \gamma) \phi_i - \delta \dot{\phi}_i $$

An active term is calculated from the difference of the left and right muscle activations $M_i^l$ and $M_i^r$ multiplied by a gain $\alpha$. A stiffness term is calculated from the muscle activities, the tonic

<figure><img src="image_5.png" /><figcaption>**FIGURE 2** | The CPG model. (A) The axial (i.e., trunk and tail) spinal network model is constituted by a double chain of 50 oscillators, i.e., 25 segments of which 19 are shown (green). Four oscillators (blue) control the limbs. Thicker arrows denote stronger couplings. For simulations of the isolated CPG, a randomly fluctuating tonic drive mimicking the pharmacological activation used in Ryczko et al. (2015) was applied to all oscillators. (B) Intersegmental phase lags from 10,000 simulations of the isolated CPG using different random seeds (intersegmental phase lag calculated by taking the average of intersegmental phase lags between segments 8–12, see Methods). Simulations are ordered by decreasing intersegmental phase lag on the vertical axis. A positive phase lag corresponds to a rostrocaudal traveling wave (i.e., from head to tail), a zero phase lag to a standing wave, and a negative phase lag to a caudorostral traveling wave. (C) Trimodal distribution of intersegmental phase lags. (D) Cycle durations vs. intersegmental phase lags. A linear fit was applied to the dataset. The square of the correlation coefficient and the significance of the fit are given.</figcaption></figure>
<table><caption>**TABLE 2** | CPG coupling parameters.</caption><thead><tr><th>Coupling type</th><th>Strength w<sub>ij</sub></th><th>Phase bias φ<sub>ij</sub> (rad.)</th></tr></thead><tbody><tr><td>Intersegmental, rostrocaudal</td><td>5</td><td>0.066:2π</td></tr><tr><td>Intersegmental, caudorostral</td><td>1</td><td>-0.066:2π</td></tr><tr><td>Intrasegmental, lateral</td><td>10</td><td>π</td></tr><tr><td>Interlimb, rostrocaudal</td><td>3</td><td>π</td></tr><tr><td>Interlimb, caudorostral</td><td>30</td><td>π</td></tr><tr><td>Interlimb, lateral</td><td>10</td><td>π</td></tr><tr><td>Limb to axial oscillators</td><td>30</td><td>4 (BTS: 5.5*)</td></tr><tr><td>Axial to limb oscillators</td><td>2.5</td><td>-4 (BTS: -5.5*)</td></tr></tbody></table>

<table><caption>**TABLE 3** | Other CPG parameters.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value (mean ± SD)</th></tr></thead><tbody><tr><td>Amplitude convergence factor</td><td>a</td><td>5</td></tr><tr><td>Saturation rate</td><td>b</td><td>500</td></tr><tr><td>Drive random walk convergence factor</td><td>c</td><td>0.001 (in vitro)<br>0 (in vivo)</td></tr><tr><td>Drive random walk step size</td><td>σ</td><td>0.03 (in vitro)<br>0 (in vivo)</td></tr><tr><td>Drive</td><td>d<sub>i</sub></td><td>See Table 4</td></tr><tr><td>Saturation threshold</td><td>d<sub>i</sub><sup>th</sup></td><td>0.3 (axis, in vitro)<br>3 (axis, in vivo)<br>0.09 ± 0.02 (limbs, in vitro)<br>1.27 ± 0.02 (limbs, in vivo)</td></tr><tr><td>Excitability</td><td>e<sub>i</sub></td><td>1.1 ± 0.07 (axis)<br>0.8 ± 0.05 (forelimbs)<br>0.5 ± 0.03 (hindlimbs)</td></tr><tr><td>Axial proprioceptive feedback, ipsilateral</td><td>w<sup>ipsi</sup></td><td>See Table 4</td></tr><tr><td>Axial proprioceptive feedback, contralateral</td><td>w<sup>contra</sup></td><td>-w<sup>ipsi</sup></td></tr></tbody></table>

stiffness $\gamma$, a stiffness gain $\beta$ and the joint angle $\phi_i$. A damping term is calculated from a damping constant $\delta$ and the joint angular velocity $\dot{\phi}_i$. Parameter values are given in **Table 5**.

In simulations, a delay of 10 ms was introduced between the CPG outputs $x_i$, $x_{i+25}$ and the corresponding muscle activations $M_i^l, M_i^r$, respectively, as a minimum to account for the muscle activation dynamics. This delay was not necessary in robot experiments since the motor torque controller already introduces a larger delay of the order of 50 ms, which is consistent with the reported range (50 ms to 1 s) of the low-pass filter properties of muscle contraction (Partridge, 1965).

## Limb Joints

For the limbs, the oscillator phase $\theta_i$ is used directly as a representation of the desired position, with a piece-wise linear transfer function that modulates the swing and stance rotation speeds such as to obtain a duty factor of 77% (Ashley-Ross and Lauder, 1997; Ashley-Ross et al., 2009). For backward terrestrial

stepping the direction of limb rotation was inverted by using $-\theta_i$ instead of $\theta_i$.

## Sensory Feedback

Proprioceptive feedback signals $s_i$ are derived from the joint angles $\phi_i$ by simulating the activity of stretch receptors: $s_i = w^{ipsi}s_i^{ipsi} + w^{contra}s_i^{contra}$, with $s_i^{ipsi}$ and $s_i^{contra}$ the positive part of $\phi_i$ and $-\phi_i$, respectively, for the left side ($-\phi_i$ and $\phi_i$ for the right side), and $w^{ipsi}$ and $w^{contra}$ the feedback weights. Since the axial part of the CPG model has 25 segments (each containing 2 hemisegmental oscillators) and the robot only 8 active axial joints (**Figure 1B**), some mapping is necessary. The signal from each joint is sent to the 3 neighboring segments, while only the middle segment is used to drive the joint muscles (**Figure 1B**). This leaves segments 3 and 16 without feedback, which is reasonable since the amplitude of the body curvature is smallest at these positions in the animal (Karakasiliotis et al., 2013).

In some simulations, an additional term was added to the $\dot{\theta}_i$ equation for limb oscillators to represent excitatory proprioceptive feedback from the limbs, as used in a previous study (Harischandra et al., 2011). Here a simplified form was used:

$$ \dot{\theta}_i = 2\pi\nu_i + \sum_j r_j w_{ij} \sin (\theta_j - \theta_i - \phi_{ij}) + w^{limb} \max \left( 0, 1 - \frac{|\phi_i - \phi_i^0|}{\frac{\pi}{2}} \right) $$

Here $w^{limb}$ is the feedback weight, $\phi_i$ the joint angle of the robot rotational limb and $\phi_i^0$ the angle at the transition from stance to swing. The feedback is maximal at the end of the stance and decreases linearly on either side until it reaches zero. The rate of decrease is such that the feedback is non-zero for half of the leg rotation. The value is always positive or zero, so this feedback term can only have an accelerating effect.

## Mechanical Simulation

3D simulations of the robot were performed using the Webots 6 software (Cyberbotics, Switzerland), which is based on the Open Dynamics Engine (ODE, www.ode.org). The physics engine was extended with a hydrodynamics model that includes reactive and resistive forces (Porez et al., 2014). The passive tail fin was modeled as a chain of 10 small segments with passive stiffness. The physics was simulated with a time step of 0.5 ms. The robot controller used a time step of 1 ms, the minimum value supported by Webots. This was just too high for a stable simulation of the muscles, so the physics plugin was used to implement the muscle model and set the joint torques.

## Robot Hardware

The robot *Salamandra robotica II* (Crespi et al., 2013) is made of a head module (9.6 cm long), 8 active modules (9.6 cm long each) and a 24.6 cm long, passive, flexible caudal fin (**Figure 1A**). This allowed the robot to approximately reach the tail length/total body length ratio of the real animal (around 0.5–0.6, see Ryczko et al., 2015). Each module actuates an axial joint with motion

<footnote>*BTS, backward terrestrial stepping.</footnote>
<table><caption>**TABLE 4** | Parameters regulating the CPG activity.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Drives d<sub>i</sub></th><th colspan="2">Feedback weights</th></tr><tr><th>Seg. 1-3</th><th>Seg. 4-25</th><th>Limbs</th><th>w<sub>psi</sub></th><th>w<sub>limb</sub></th></tr></thead><tbody><tr><td colspan="6"><strong>Robot experiments with 5 individuals (Figure 5)</strong></td></tr><tr><td>Swimming</td><td></td><td>1</td><td></td><td>10</td><td>0</td></tr><tr><td>Forward terrestrial stepping</td><td>0.60 ± 0.02</td><td></td><td>1.00 ± 0.04</td><td>0</td><td>0</td></tr><tr><td>Forward underwater stepping</td><td>0.42 ± 0.01</td><td></td><td>0.71 ± 0.03</td><td>0</td><td>0</td></tr><tr><td>Backward terrestrial stepping</td><td>0.23 ± 0.01</td><td></td><td>0.46 ± 0.02</td><td>0</td><td>0</td></tr><tr><td>Struggling</td><td>0.27 ± 0.01</td><td></td><td>0.38 ± 0.01</td><td>0</td><td>0</td></tr><tr><td colspan="6"><strong>Other robot experiments</strong></td></tr><tr><td>Swimming without regulation (Supplementary Movie 1)</td><td></td><td>1</td><td></td><td>0</td><td>0</td></tr><tr><td>Swimming with differential drive (Supplementary Movie 2)</td><td>0.9</td><td></td><td>1</td><td>0</td><td>0</td></tr><tr><td>Forward terrestrial stepping with feedback (Figure 6E)</td><td></td><td>0.61</td><td>0.98</td><td>from<br/>-10 to 6</td><td>0</td></tr><tr><td colspan="6"><strong>Simulations</strong></td></tr><tr><td>Isolated CPG (simulation of in vitro experiments) (Figures 2, 3)</td><td></td><td>0.1 ± 0.01</td><td></td><td>0</td><td>0</td></tr><tr><td>Swimming without regulation (Supplementary Figure 4A)</td><td></td><td>1.34</td><td></td><td>0</td><td>0</td></tr><tr><td>Swimming with differential drive (Supplementary Figure 4B)</td><td>1.03</td><td></td><td>1.34</td><td>0</td><td>0</td></tr><tr><td>Swimming with axial feedback (Supplementary Figures 4C,D)</td><td></td><td>1.34</td><td></td><td>21</td><td>0</td></tr><tr><td>Forward terrestrial stepping without regulation (Figure 6A)</td><td></td><td>0.98</td><td></td><td>0</td><td>0</td></tr><tr><td>Forward terrestrial stepping with differential drive (Figure 6B)</td><td></td><td>0.63</td><td>0.98</td><td>0</td><td>0</td></tr><tr><td>Forward terrestrial stepping with axial feedback (Figure 6C)</td><td></td><td>0.98</td><td></td><td>-0.65</td><td>0</td></tr><tr><td>Forward terrestrial stepping with limb feedback (Figure 6D)</td><td></td><td>0.98</td><td></td><td>0</td><td>3.7</td></tr></tbody></table>

Robot experiments used lower drives for swimming to stay in the robot operating range (see Results). Standard deviations for the drive reflect variations between the simulated individuals. Drive values are shown centered across two or three columns in cases where the same value was applied to the corresponding groups of oscillators.

<table><caption>**TABLE 5** | Muscle parameters.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Muscle active gain (N-m)</td><td>α</td><td>0.4 (BTS, ST: 4*) (simulation)</td></tr><tr><td></td><td></td><td>0.5 (BTS, ST: 5*) (robot)</td></tr><tr><td>Muscle stiffness gain (N-m/rad)</td><td>β</td><td>1.2 (BTS, ST: 12*)</td></tr><tr><td>Muscle tonic stiffness (no unit)</td><td>γ</td><td>0.2</td></tr><tr><td>Muscle damping (N-m-s/rad)</td><td>δ</td><td>0.1</td></tr></tbody></table>

restricted to the horizontal plane; the two girdle modules also include rotational joints for the limbs. The entire robot measures 111 cm and weighs 2.48 kg. The robot modules have LEDs on the covers which were used to track the robot's motion with two Basler A622F video cameras (15 frames/s) to cover the whole track length (6 m) with an accuracy of ± 1 cm. The Supplemental Movies of the robot were captured with another camera at 15 or 30 frames/s. Two adaptations were made to the robot to reproduce the different behaviors. During forward underwater stepping, the buoyancy was adjusted by adding 72 g of lead in the head. This corresponds to +2.9% of the total robot weight, or + 41.8% of the normal weight of the head module (172 g). During struggling, tape was added under the feet to increase slipping, mirroring the conditions of the animal experiments (Ryczko et al., 2015).

## Distributed Electronics and Control Software

The robot controller is distributed: each module reads the position and velocity of the local joint and computes the control loop for the corresponding part of the CPG and the joint's virtual muscles, with a time step of 10 ms. The numerical integration of CPG segments with floating-point operations required a modification to the hardware described by Crespi et al. (2013): the modules were upgraded to include an LPC2129 ARM7TDMI microcontroller running at 60 MHz, as already present in the head. Communications between modules are restricted to drive signals from the head and CPG couplings between adjacent modules, sent over the CAN bus running at 1 Mbps. The leg positions are set by PD controllers using the motor encoders. The axial torques are set by PI controllers using current sensing. The CPG state from each module was recorded by logging coupling and debug messages sent over the CAN bus. This logging was done on an external computer, by spying on the bus using long, thin wires attached to the caudal end of the robot. Two modifications were made to the distributed controller between the initial tests and the final version (see Results): (1) The numerical integration of the CPG was changed to estimate the phases $θ_j$ and amplitudes $r_j$ of coupling sources at the time of integration using a linear extrapolation of the values from the two latest CAN messages and their times of arrival (coupling terms are dropped entirely from the integration if the two previous CAN messages are older than 200 or 400 ms, respectively); (2) The phases and amplitudes were encoded in

<footnote>*BTS, backward terrestrial stepping; ST, struggling.</footnote>
CAN messages as 16-bit half-floats rather than 32-bit floats, such that coupling signals from a module to a particular neighbor would fit in a single message.

## Selection of Parameter Values

The CPG parameters were hand-tuned to reproduce the distribution of intersegmental phase lags and cycle durations from *in vitro* recordings of isolated spinal cords (Ryczko et al., 2015). To emulate the motor behaviors displayed by the animal *in vivo*, the CPG model was subjected to higher excitatory drives which were tuned to reproduce typical electromyographic patterns for each of the five motor behaviors, in terms of cycle frequency and intersegmental phase lags (Ryczko et al., 2015). However, the cycle frequencies targeted with the model were set to half that of the *in vivo* recordings. This was chosen to reflect the scaling of locomotion frequency with body mass observed in animals (Bejan and Marden, 2006). For robot experiments, the target frequency for swimming was further lowered to 1.1 Hz due to the limits of operation of the robot (in particular torque limits). Other than the drive levels and sensory feedback, the only changes from *in vitro* to *in vivo* CPG conditions were in the average saturation thresholds which had to be increased to match the higher drives used *in vivo*.

The parameter space for the virtual muscles and proprioceptive feedback was explored systematically using the 3D mechanical simulation of the robot during swimming. An “average individual” was used by setting the standard deviations of the CPG excitabilities to zero, to increase reproducibility (this restriction was relaxed for robot experiments). We used uniform muscle parameter values for trunk joints, and progressively smaller values in the tail to emulate body taper: the values of $\alpha$ and $\beta$ in modules 6, 7, and 8 were multiplied by a factor 0.7, 0.5, and 0.2, respectively. The same feedback parameter values were used for all joints, and the same feedback weights (with opposite signs) were used for ipsilateral and contralateral projections. Initial tests were made with a tonic muscle stiffness $\gamma = 0$: this parameter is mostly redundant with the stiffness gain $\beta$ for a given (non-zero) amplitude of CPG oscillations. A uniform excitatory drive was used for all oscillators, which in absence of feedback results in high intersegmental phase lags inappropriate for swimming. The drive was set to 1.34, corresponding to a swimming frequency of 1.47 Hz (in absence of sensory feedback), which is close to our target of half the frequency observed in the animal (2.78 Hz and 3.12 Hz during EMG and kinematic recordings, respectively, Ryczko et al., 2015). We selected muscle and feedback parameter values that showed a significant increase in swimming speed and high stability of the CPG and kinematic patterns, while keeping the joint torques close to the robot’s limit of 0.7 Nm (Supplementary Figures 1–3).

The five salamander behaviors were first reproduced in simulation without sensory feedback, by tuning the CPG drive levels and optimizing the limb-body phase bias for speed of locomotion. Additional simulations were done with varying strengths of axial proprioceptive feedback during forward terrestrial stepping and swimming. Further simulations were made with proprioceptive limb feedback during forward terrestrial stepping.

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

**FIGURE 3** | The CPG spontaneously switches between axial motor patterns as a function of a fluctuating background drive strength. The strength of the simulated pharmacological drive applied to the CPG is shown together with the outputs of the oscillators of the left forelimb, left trunk segments 8 and 10, and the left hindlimb. Before the switch, limb oscillators (blue lines) are saturated by the drive strength. Therefore, they show low amplitude oscillations and are entrained to the higher frequency of the trunk oscillators, and the motor pattern in the trunk segments (green) is a rostrocaudal wave (white dots). Then, a progressive decrease in drive strength occurs through random fluctuations, and this progressively de-saturates the limb oscillators. The de-saturation allows limb circuits to oscillate at higher amplitude, causing a switch (arrows) from a fast rostrocaudal wave to a slower caudorostral wave (black dots). After the switch, the de-saturated limbs show high amplitude oscillations and therefore entrain to their slower frequency the trunk oscillators, because effective connection strength in the model is proportional to the amplitude of oscillations from the sender (see Methods). The same switches have been observed in the isolated spinal cord [Figures 5A,B of Ryczko et al. (2015)].

The limb-body phase bias was optimized again on the robot, due to the different zero-point reference and the backlash in the gears. The main simulation results were then reproduced on the robot with five different “individuals,” which were modeled by initializing the CPG parameters using different random seeds. The descending drives were adjusted for each “individual” and each motor behavior. Other parameters were sometimes adjusted between behaviors but always using the same values for all individuals (see Results). The movies shown in the Supplemental Materials were prepared using an average individual.

## Data Processing

The joint angles from simulations and robot experiments were calculated in Matlab by fitting the kinematic chain of the robot to the positions of the LEDs. The CPG and kinematic intersegmental phase lags were calculated in Matlab from the CPG output and joint angle oscillations using the same algorithm: The timing of each cycle was defined as the centroid of the positive part of each oscillation (Knüsel et al., 2013). These timings were used to calculate a median phase lag (over time) for each pair of consecutive segments (CPG lag) or consecutive joints (kinematic lag). The CPG intersegmental phase lag was calculated using the average of the median phase lags between segments 8 to
12. For simulations, the kinematic intersegmental phase lag was calculated using the average of the median phase lags between joints 3 to 5. For robot experiments, joint 5 was often an outlier due to the robot torque limits, so joints 2 to 4 were used instead.

## Statistics

Data are given as means ± standard deviation (SD) unless specified otherwise. Correlations between variables were evaluated in SigmaPlot 11.0 using the Pearson Product Moment Correlation test.

## RESULTS

### The Isolated CPG Model Reproduces the Main Features of Recordings From Isolated Spinal Cords

We found that asymmetric intersegmental couplings, together with different excitabilities and saturation thresholds between forelimb, hindlimb and axial oscillators, enabled the spinal cord model to reproduce the trimodal distribution of phase lags observed *in vitro* during fictive locomotion (Ryczko et al., 2015). The values shown in Tables 2, 3 were found to produce a range of phase lags similar to the biological data, with peaks centered on 6.7 ± 1.3, 1.3 ± 1.3, and -6.1 ± 1.0% (Figures 2B,C). The positive correlation between phase lag and cycle duration (Ryczko et al., 2015) was also reproduced (Figure 2D). Furthermore, the small fluctuations in the excitatory drive over time allowed the model to reproduce the spontaneous switches between slow caudorostral waves and fast rostrocaudal waves of axial activity reported in the isolated spinal cord (Delvolvé et al., 1999; Ryczko et al., 2015). Figure 3 shows an example of the model producing such a switch.

### Improving Swimming With Proprioceptive Feedback Requires Specific Muscle Stiffness and Damping Properties

Using a uniform drive yielded an intersegmental phase lag of 6.6% in an average individual (leftmost peak of the distribution). This resulted in inefficient swimming, with too many nodes in the traveling wave (Supplementary Movie 1). We looked for muscle parameters that would allow proprioceptive feedback to improve swimming by decreasing the phase lag toward more physiological values. For the active gain $α$, a value of 0.4 proved optimal, as higher values (together with higher stiffness $β$ or $γ$) would have given higher swimming speeds but would have required torques beyond our robot's limits. Systematic tests in the $β, δ, w^{ipsi}$ space showed a single region where feedback increased the swimming speed thanks to a decrease of the phase lag, while keeping the CPG rhythm stable. This stable region corresponds to a stiffness gain $β$ between 1.6 and 2.3, a damping $δ$ between 0.05 and 0.15 and feedback weights $w^{ipsi} = -w^{contra}$ between 17 and 22 (Supplementary Figures 1-3). Further tests in the $β, γ, w^{ipsi}$ space showed that we could trade some fitness gain $β$ for tonic stiffness $γ$. We settled on $β = 1.2, γ = 0.2$, which give qualitatively reasonable passive mechanical properties (Table 5).

### Simulation Results Transferred to the Robot Following Some Adaptations

The simulation results could be reproduced on the robot, with the following changes made based on qualitative judgements: the muscle active gain $α$ had to be increased from 0.4 to 0.5 to obtain reasonable amplitudes of oscillation during swimming. Uniform muscle parameter values were used in all robot joints: the tapering of the active and stiffness gains was removed to obtain reasonable amplitudes of oscillations in the tail and good swimming speeds. Due to the limits of operation of the robot, the target frequency for swimming had to be lowered down to 1.1 Hz<fnref n="2" />. The feedback weights were reduced to 10. Sensory feedback from the neck joint was removed as it was destabilizing, leading to aperiodic rhythms. Finally, in simulation we found that a common limb-body phase bias gave near-optimal speed for all stepping behaviors. This was not the case with the robot: a specific limb-body phase bias was required for backward stepping (Table 2) to obtain the optimal speed for that behavior.

### Two Drive Signals Suffice to Reproduce the Five Motor Behaviors, but Backward Terrestrial Stepping and Struggling Require Stronger Muscle Contractions

We found that tonic drives with only two different values applied to different parts of the CPG were sufficient to reproduce qualitatively the five motor behaviors with the robot, as shown by movies (Supplementary Movies 2-6), frame sequences (Figure 4) and robot kinematics (Figure 5, Supplementary Figure 4B). In particular, the differences in CPG and kinematic intersegmental phase lags between the five behaviors were reproduced (Figures 5F,G), as well as the differences in cycle durations (Figure 5H). Figures 5A-E show the CPG outputs and joint oscillations for a single individual. The rostrocaudally increasing lag between CPG and kinematic waves is reproduced (increasing gap between the thick red and thick black lines).

Swimming could be obtained by sending a strong drive (i.e., saturating limb oscillators, Hypothesis 3) to the whole CPG, with a slightly lower drive to the most rostral oscillators (segments 1-3 in Figure 2) to adjust the phase lag as proposed in numerical simulation of the lamprey locomotor CPG (Kozlov et al., 2009). The other behaviors were obtained by adjusting the drive to the limb oscillators independently from the drive to the axial oscillators (Figure 5H, Supplementary Movies 3-6). Higher axial phase lags required a greater relative difference between the two drives, and higher frequencies required higher values of both drives.

While two drives were sufficient to generate the CPG activity patterns for all behaviors, we found that stronger forces from the virtual muscles were required to reproduce the

<footnote n="2">²This meant producing a swimming gait in the frequency range of forward terrestrial stepping, i.e. with drives under the limb saturation threshold, so limb oscillators had to be silenced artificially.</footnote>
<figure><img src="image_7.png" /></figure>

**FIGURE 4** | Frame sequences showing from top to bottom one complete stride for the five motor behaviors. Light-emitting diodes were tracked on each module for kinematic analysis. A black dot was added on each frame to illustrate robot progression. Inter-frame time intervals in ms are, respectively: 200 (swimming, SW); 133.3 (forward underwater stepping, FUS); 166.7 (forward terrestrial stepping, FTS); 400 (backward terrestrial stepping, BTS); 466.7 (struggling, ST). Scale bar (white), 10 cm. The different background colors are due to the different environments: water tank for SW and FUS, wooden board for the other behaviors. The behaviors are close to those observed in the real animal [Figure 1 of Ryczko et al. (2015)].

kinematics of backward terrestrial stepping and struggling. To avoid introducing additional parameters, this was implemented by increasing the muscle gains α and β. A 10-fold increase was found appropriate to avoid large deviations between the CPG activity and the kinematics (Table 5, Figures 5A,E).

## Proprioceptive Feedback Can Regulate the Phase Lag During Swimming and Forward Terrestrial Stepping, and Reduces Variability

With axial proprioceptive feedback, swimming could be reproduced with a single drive to the whole CPG (Figures 5B,I,
**Supplementary Movie 8).** In absence of feedback, the axial network of the five robot "individuals" produces intersegmental phase lags of 6.1 ± 1.4% (Supplementary Figure 4A, Supplementary Movie 1). With feedback, this could be reduced to 1.9 ± 1.1%, which matches the values observed in the animal (1.89 ± 0.25%, Ryczko et al., 2015). Feedback also reduced the variability between individuals: Without feedback, the individual corresponding to the leftmost curve in Figure 5I (dashed lines) was an outlier. Using feedback with identical weights in all individuals (solid lines), the outlier showed phase lags similar to the other individuals.

Axial proprioceptive feedback could also replace differential drives (i.e., different drive values for different parts of the CPG) as a means of obtaining axial phase lags close to zero in simulations of forward terrestrial stepping. However, this required feedback weights of opposite signs compared to swimming. Using uniform drives and no feedback, the CPG generated negative phase lags (Figure 6A). With uniform drives, positive values of $w^{ipsi}$ (as used in swimming) had a destabilizing effect on the CPG, with oscillators failing to lock in frequency. With differential drives as in Figure 6B, we found that positive values of $w^{ipsi}$ were counterproductive, decreasing again the CPG phase lag to negative values. However, negative values of $w^{ipsi}$ produced higher phase lags as desired and could be used to reproduce forward terrestrial stepping with uniform drives (Figure 6C). This effect of an increasing intersegmental phase lag with decreasing feedback weights during forward terrestrial stepping was reproduced in robot experiments (Figure 6E).

Instead of axial feedback, the differential drive could also be replaced with limb proprioceptive signals fed back to the limb oscillators. This was only tested in simulation (Figure 6D). The effect of feedback here was again to increase the negative CPG phase lags toward slightly positive values, close to zero.

## A Passive Tail Decreases the Drag During Underwater Stepping

During forward underwater stepping, the locomotor performance could be improved by using “passive” tail segments (Supplementary Movie 7), similar to the animal which shows passive tail undulations during this locomotor behavior (see Cabelguen et al., 2014). This was implemented using a third drive level of zero to the tail segments. The tail CPG was then inactive, but the robot modules continued to generate torques corresponding to the passive parts of the muscle model (the terms that remain when $M_i^l$ and $M_i^r$ are 0). Measurements of the hydrodynamic forces in the 3D mechanical simulation suggest that for this particular gait, the undulations of passive tail segments allow the caudal fin to generate more thrust than in the case of active tail segments (Figure 7). The drag at the head and girdles is also reduced in the passive case. The net drag on the body axis is thus reduced from -0.0422 ± 0.0003 N to -0.0295 ± 0.0004 N (standard errors).

## Coupling Delays Introduce Systematic Phase Biases in the Distributed Robot Controller

Initial tests on the robot with the distributed controller gave non-uniform phase lags along the body, unlike what was seen in simulation. We investigated the issue using a chain of 7 simple modules (no girdles) and a CPG model with symmetrical ascending and descending coupling weights and phase biases of 5%. We found increased phase lags between the first modules and decreased phase lags between the more caudal ones (Figure 8A). An analysis of the coupling terms used in the numerical integration of the CPG showed that rostral modules were significantly slowed down by caudal ones (Figure 8B). This suggested that the $\theta_j$ values (the phases of the couplings' source oscillators) used in the numerical integration of the target oscillators were out of date. We modified the numerical integration to estimate the state of the source oscillator at the time of integration using a linear extrapolation of the two previous coupling messages and their times of arrival (see Methods). This considerably reduced the slow-down effect and yielded almost uniform phase lags (Figures 8C,D). Further improvements (not shown) were obtained by encoding the coupling phases and amplitudes as 16-bit half-floats. This halved the number of messages sent over the CAN bus and helped decrease the rate of transmission errors.

## The Distributed Controller Allows for Autonomously Moving Robot Parts

We found that the distributed implementations of the CPG and muscle models have the side effect of making the robot modular at runtime. We conducted forward terrestrial stepping experiments with screws between some modules removed, causing the robot to break into parts (no other changes were made to the hardware or software). Each part kept functioning, still coordinated by its own section of the CPG model. This is demonstrated in Supplementary Movie 9, Figure 9.

## DISCUSSION

### A Modulable CPG Architecture

Our results suggest that the answer to question 1 is yes: a modulable CPG provides a robust framework for generating multiple motor patterns, such that different motor behaviors do not necessarily require dedicated CPGs. This concept was proposed by Grillner (1981) as the “unit burst generator” theory, which states that independent rhythmogenic circuits can be flexibly coupled from one behavior to another. Such circuits have been identified in many animals. In the salamander, specific spinal hemisegments have been shown to control muscles of the trunk (Ryczko et al., 2010), tail (Charrier and Cabelguen, 2013) and limbs (Cheng et al., 1998; Lavrov and Cheng, 2004; Ijspeert et al., 2007). Other examples include the spinal hemisegments in the lamprey (Cangiano and Grillner, 2003, 2005; Cangiano et al., 2012), the flexor
<figure><img src="image_8.png" /></figure>

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

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

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

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

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

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

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

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

**FIGURE 5** | CPG output and joint kinematics during the five motor behaviors reproduced with the robot. Results obtained using a single drive level and proprioceptive feedback for swimming, two drive levels (to limb and axial oscillators) and no feedback for the other behaviors, with a third drive of zero for the tail during underwater stepping. (A–E) and (F,G) emulate the biological data illustrated in Figures 2, 3 of Ryczko et al. (2015), respectively. (A–E) Kinematic angular oscillations (thin black (Continued))
**FIGURE 5** | lines) and CPG outputs $x_i$ (thin red lines) are shown for each joint, with circular markers indicating the centroid of the positive half of each cycle. Rostrocaudal (A-C, line descending to the right), standing (D, almost vertical line) and caudorostral (E, line descending leftward) kinematic waves (thick black lines) followed the CPG activity (thick red lines) with variable delays during the 5 motor behaviors: backward terrestrial stepping (BTS), forward swimming (SW), forward underwater stepping (FUS), forward terrestrial stepping (FTS) and struggling (ST). All figures are from the same simulated individual. (F,G) Intersegmental phase lag distributions for CPG waves (F) and for kinematic waves (G) in the robot trunk. Each marker represents a single recording. (H) Comparison of the tonic drive signals applied to the CPG and resulting oscillation frequencies (Hz) for each motor behavior. Data represent mean ± SD over the 5 simulated individuals. The overlap in limb drive values between SW and FTS is due to reduced drive and limb saturation thresholds during SW to accommodate for the robot's torque limitations (See Methods). Without these limitations, the limb and axial drives for SW would be 1.34 as in simulation. (I) Rostrocaudal axial waves generated during SW by the CPG without sensory feedback (dashed lines) and with sensory feedback (solid lines), with uniform drive $d_i$ = 1, in the 5 simulated individuals. For each individual, the wave with feedback is horizontally positioned in the figure to connect to the wave without feedback. The experiments with sensory feedback correspond to the SW data in (F-H).

and extensor networks in the left and right side of the mouse spinal cord (Hägglund et al., 2013), the crayfish swimmeret system (Mulloney and Smarandache-Wellmann, 2012) and the networks controlling individual leg joints in the stick insect (Büschges et al., 1995). A modeling study of insect locomotion suggests that the recruitment of a single neural structure for various behaviors also applies to situations where locomotion is largely driven by sensory feedback (Schilling and Cruse, 2020).

## Oscillator Couplings

Biological data indicate that in salamanders, limb activity can occur together with traveling waves *in vivo* and *in vitro* (Ryczko et al., 2015). In our previous robotic study (Ijspeert et al., 2007), limbs projected to all axial oscillators (forelimb oscillators to trunk oscillators, and hindlimb oscillators to tail oscillators). Two axial outputs were therefore possible: either a standing wave when limb oscillators were active (during stepping), or a rostrocaudal wave when limbs were saturated (during swimming). Here we made the limb oscillators project only to neighboring axial oscillators (Hypothesis 1), which gives more flexibility for the coordination of axial oscillators when limbs are rhythmically active (Ijspeert et al., 2005; Knüsel et al., 2013). In the animal, a unidirectional connection from excitatory interneurons generating the limb rhythm to those generating the axial rhythm would be sufficient to impose the slow limb oscillations to the axial segment, according to a detailed model of salamander spinal networks based on Hodgkin-Huxley neurons (Bicanski et al., 2013).

In our model, the asymmetry between ascending and descending coupling weights $w_{ij}$ is required to reproduce a wide diversity of axial phase lags with consistent values along the axis as observed *in vitro* and *in vivo* (Ryczko et al., 2015). Symmetric weights as used in Ijspeert et al. (2007) produce non-uniform phase lags along the axis when the oscillators have different intrinsic frequencies. Future studies should determine whether such a coupling is present in salamanders and how it is implemented. Possibilities include an asymmetry at the neuroanatomical level (dominance of descending projections, see Buchanan et al., 1989; Buchanan, 2001 in lamprey) or in electrophysiological terms (stronger synaptic strengths toward caudal segments, see Smarandache et al., 2009 in crayfish; more spikes per locomotor cycle in neurons projecting caudally, see Mulloney et al., 2006 in crayfish).

## Oscillator Frequencies and Saturation

For our CPG model to be able to generate the three types of axial waves recorded *in vitro* and *in vivo* in salamanders, and the positive correlation between cycle duration and phase lag (Ryczko et al., 2015), we had to modify the intrinsic frequency of limb networks compared to our previous study (Ijspeert et al., 2007). Forelimb and hindlimb oscillators still have an intrinsic frequency slower than axial oscillators, but here forelimb oscillators are faster than hindlimb ones (Hypothesis 4). Data in mammals suggest that forelimbs deprived of normal interactions with the hindlimb networks tend to accelerate *in vivo*. Indeed, in adult cats where the spinal cord is partially lesioned, forelimb and hindlimb rhythms often dissociate, and forelimbs adopt a faster rhythm, yielding a 2:1 forelimb-hindlimb coupling (for review, see Frigon, 2017). At the cellular level, modification of a single conductance controlling burst termination should be sufficient to make limb segments generate slower oscillations, as suggested by a detailed Hodgkin-Huxley model of a salamander spinal segment (Bicanski et al., 2013).

A hypothesis that we kept from our previous work (Ijspeert et al., 2007) is that with a strong descending drive, limb networks "saturate" whereas axial oscillators do not (Hypothesis 3). Future studies should examine whether and how such a function is implemented in the animal. It could be a differential recruitment of specialized interneuron populations as a function of drive strength, as documented as a function of speed in zebrafish (McLean et al., 2007, 2008; Gabriel et al., 2011; Ampatzis et al., 2014, for review see Berg et al., 2018) and mice (Talpalar et al., 2013, for review see Kiehn, 2016). It could also involve a shift in the active set of reticulospinal neurons as a function of speed/gait. Some reticular neurons increase their firing specifically during swimming in salamanders (Lowry et al., 1996). Different reticulospinal neurons are activated as a function of speed in zebrafish (Kinkhabwala et al., 2011).

## Regulation Through Descending Drives

In answer to question 2, our results (Figure 5, Supplementary Figure 4B, Supplementary Movies 2–6) suggest that independent drive levels to a few parts of a CPG network (here two, or three to reproduce passive tail undulations) are sufficient to emulate a diversity of motor behaviors. In the model, the regulation of CPG activity by descending drives can be understood intuitively. The drive signals control the intrinsic (uncoupled) frequencies of the oscillators. Because rostrocaudal couplings are stronger than caudorostral couplings,
<figure><img src="image_17.png" /><figcaption>**FIGURE 6** | Effect of axial or limb proprioceptive feedback on CPG activity and kinematics during forward terrestrial stepping (FTS). (A–D) Joint kinematics (black) and CPG activity (red, dashed for limb oscillators) in simulation. Circular markers indicate the centroid of the positive half of each cycle. The limb-body phase bias $\varphi_{ij}$ was adapted in each case for optimal speed of locomotion. (A) With uniform drive ($d_i = 0.98$) and no feedback, the CPG with active (non-saturated) limb oscillators produces caudorostral waves of activity. (B) Standing waves of CPG activity can be obtained by using a different drive of 0.63 for the limb oscillators. (C) Standing CPG waves can also be obtained using uniform drives and axial proprioceptive feedback, with $w^{\text{ipsi}} = -w^{\text{contra}} = -0.65$. (D) Standing waves of CPG activity could also be obtained with uniform drives using proprioceptive limb feedback, with $w^{\text{limb}} = 3.7$. (E) Effect of the axial feedback weight on CPG axial intersegmental phase lag in robot experiments. Before introducing feedback, differential drives to the limb ($d_i = 0.98$) and axial ($d_i = 0.61$) oscillators were used to increase the phase lag toward zero. Positive ipsilateral feedback weights (as those used during swimming with sensory feedback, **Figure 5I**) decreased the phase lag, whereas negative ipsilateral weights increased it.</figcaption></figure>

a segment will entrain a slower or faster caudal neighbor, and the resulting common frequency will be close to the frequency of the rostral segment. However, the faster segment will lead the slower one with a delay that increases with the difference in uncoupled frequencies (this delay being in addition to the coupling's natural phase bias). This effect will propagate down the chain of segments, such that the resulting frequency and phase lag of the whole chain can be controlled by adjusting two

values: the uncoupled frequency of the first segment, and that of the other segments.

This mechanism of regulation is close to the “trailing oscillator hypothesis,” which states that the oscillator of higher excitability becomes the leader of the chain. This hypothesis is based on lamprey experimental data showing that increasing the excitability in caudal segments causes a switch from a rostrocaudal to a caudorostral wave in isolated spinal cords
<figure><img src="image_18.png" /><figcaption>**FIGURE 7** | Passive tail CPG segments reduce drag during forward underwater stepping (FUS) in simulation. The hydrodynamic forces acting on different parts of the robot modules along the direction of motion are given. The caudal fin provides more thrust when the tail is passive. Forces on each module were measured for an "average individual" during 15 s after a warm-up period of 15 s, and averaged over the 30,000 time steps. The procedure was repeated 100 times with different starting conditions. The bars show the means over the 100 repetitions. Error bars for the standard error are shown, but barely visible (all standard errors are < 0.00013).</figcaption></figure>

(Matsushima and Grillner, 1990, 1992). However, this lamprey model assumed symmetrical rostrocaudal and caudorostral couplings, while we found that the coupling asymmetry is important to maintain a uniform phase lag along the chain of oscillators. A later lamprey modeling study with a detailed neural network of Hodgkin-Huxley neurons showed that dominant descending couplings allow for flexible control of forward and backward swimming with constant phase lag along the spinal cord at different speeds: the frequency and intersegmental phase lag can be controlled by adjusting the excitatory drive of the first segments compared to the remaining ones (Kozlov et al., 2009).

In our salamander model, the differential excitation of the first segments can be realized through the strong connections from the forelimb oscillators (when they are active). The regulation of the axial CPG pattern is then achieved by adjusting the excitation of the limb oscillators compared to the axis, instead of the first axial segments compared to the others as in the swimming case. This mechanism of regulation has been investigated with abstract oscillators and validated with a more detailed integrate-and-fire model (Knüsel et al., 2013).

The coordination of limb muscles was beyond the scope of this study: the limbs of our robot have a single rotational degree of freedom, and the direction of rotation was artificially inverted for backward stepping. We expect that more drives would be required in a model with more realistic limbs. Turning was also not investigated here but can in principle be obtained during swimming and stepping using different drives for axial oscillators on the left and right sides (Ijspeert et al., 2007).

## The Regulation Mechanism in the Isolated CPG

The mechanism of regulation described above, together with the differences in excitability and saturation thresholds between

forelimb, hindlimb and axial oscillators (Table 3), enable the isolated CPG model to reproduce the trimodal distribution of phase lags observed in vitro: In the model, hindlimb oscillators are intrinsically slower than forelimb oscillators. Given the random nature of the saturation thresholds, forelimbs or hindlimbs can selectively saturate due to slightly different threshold values. When all oscillators are active, the hindlimb oscillators slow down the forelimbs, and the strong local connections from limb to axial oscillators slow down the girdle segments, leading to a highly negative phase lag in the trunk and tail axial networks. This corresponds to the rightmost peak of the distribution (i.e., negative lags, Figures 2B-D). When the hindlimb oscillators saturate, the forelimb oscillators accelerate a bit but continue to slow down the first segments, yielding the phase lags that make up the middle peak of the distribution (i.e., near zero lags, Figures 2B-D). When all limb oscillators saturate, the axial network is no longer influenced by limb network activity and generates the higher, positive phase lags found in the leftmost peak (i.e., positive lags, Figures 2B-D).

This mechanism also explains the spontaneous switches between slow caudorostral waves and fast rostrocaudal waves of axial activity: In the isolated CPG model, the transitions between the active and saturated states are triggered by small fluctuations in the excitatory drive (Figure 3), which represents tonic pharmacological excitation as in Ryczko et al. (2010, 2015) or Delvolvé et al. (1999). The progressive saturation of the limb oscillators causes their oscillation amplitude to diminish as the cycle frequency increases. The model thus suggests that limb burst amplitude *in vitro* should be higher during slow caudorostral wave of activity than during a rostrocaudal wave (Figure 3).

## Regulation Through Proprioceptive Feedback

Recordings from isolated spinal cords show much more variability among salamander individuals than EMG recordings of intact animals (Ryczko et al., 2015). In response to question 3, our results from robot experiments with five swimming "individuals" suggest that local sensory feedback could explain this reduction of variability from the *in vitro* to the *in vivo* condition (Figure 5I): sensory feedback made it unnecessary to tune the drive levels in each individual (compare the standard deviation of the drives in Figure 5H for swimming vs. the other behaviors). Results from robot experiments and simulations also suggest that local sensory feedback can replace differential drive as a modulator of the CPG activity to produce forward terrestrial stepping (Figure 6) and swimming (Figures 5B,I, Supplementary Figures 4A,C, Supplementary Movies 1, 8), which answers the other part of question 3.

The regulation of our CPG model by proprioceptive feedback can be explained with the same mechanism as regulation by different drive signals. Sensory feedback has been previously reported to increase the locomotion cycle frequency through an excitatory effect on the lamprey CPG activity (e.g., Kiemel and Cohen, 2001). In our model, the addition of proprioceptive
<figure><img src="image_19.png" /></figure>

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

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

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

**FIGURE 8** | Systematic phase biases due to the slow-down effect of coupling delays in the robot's distributed controller. Undesired phase lags in the distributed controller were analyzed using a chain of 7 simple (non-girdle) robot modules and a CPG model with symmetrical ascending and descending coupling weights and a uniform intersegmental phase bias of 5%. (A) CPG intersegmental phase lags calculated from the middle segment of each module, as percentage of the cycle duration. The observed values are higher than the target of 5% in the rostral modules and lower than the target in the caudal ones. (B) Average effect of the ascending and descending couplings on the left oscillator of the middle segment in each module. Values shown correspond to the coupling terms $r_i w_{ij} \sin(\theta_i - \theta_{ij})$, averaged over the whole recording. The negative red bars show that the net effect of both coupling types is to slow down the oscillation. The effect is stronger in more rostral modules. (C) Phase lags observed after the implementation of coupling message extrapolation. The observed values are almost uniform and close to the target of 5%. (D) With coupling message extrapolation, the slow-down effect has almost vanished.

feedback in the axis increases the uncoupled frequency of the segments in the axial network. If the first segment receives no feedback, as is the case in robot experiments, its uncoupled frequency is comparatively reduced. This leads as expected to a decrease in phase lags during swimming (Figure 5I).

Interestingly, simulations showed that feedback can also regulate swimming when neck feedback is included (Supplementary Figures 1–3, 4C). This suggests that feedback has a weaker accelerating effect in the first segment than in the others, even though the feedback amplitude is comparable (see Supplementary Figure 4C). This can be explained by looking at the model equations: axial feedback adds the term $-\frac{\dot{s}_i}{r_i} \sin \theta_i$ to the instantaneous frequency $\dot{\theta}_i$ of the oscillator. The average value of this term is highly dependent on the

phase relationship between $\theta_i$ and the phase of the feedback signal $s_i$. In particular, if the kinematics follow closely the CPG output, and if we approximate $s_i$ with a sine wave, then $s_i$ will be proportional to $\cos \theta_i$. Assuming a constant amplitude $r_i$, the effect of feedback on $\dot{\theta}_i$ can be written as $k \cos \theta_i \sin \theta_i$, which averages to zero over a $2\pi$ interval for $\theta_i$. We conclude that if $\theta_i$ increases approximately linearly with time, the effect of feedback on the frequency will approach zero when the CPG-mechanical phase lag approaches zero. And this lag (the distance between the red and black dots) is indeed very small for the neck joint in Supplementary Figure 4C.

We can also explain the need for reversed axial feedback weights during forward terrestrial stepping: Excitatory axial feedback (as in swimming) accelerates the mid-trunk oscillators, which tends to decrease the intersegmental phase lag. This is
<figure><img src="image_23.png" /><figcaption>**FIGURE 9** | Frame sequences showing the behavior of the robot when split into several parts. (A) After losing its tail, the robot keeps moving (although with a malfunction in the pelvic girdle module). The tail modules continue to oscillate. The tail CPG maintains the coordination in the form of a rostrocaudal traveling wave. (B) The robot splits behind both girdles. The three robot parts keep moving independently.</figcaption></figure>

counter-productive since the unregulated intersegmental phase lag is already too low (Figure 6A). With inverted weights, axial feedback slows down the mid-trunk and increases the phase lag as desired. The inclusion of neck feedback has little importance in this case: when the limbs are active, the activity of segment 3 is largely determined by the strong connections from the forelimb oscillators, irrespective of the activity in the first and second segments. This also means that we would expect similar results with a model that includes head stabilization as observed in the animal during forward terrestrial stepping (Ryczko et al., 2015). Our results with inverted weights are reminiscent of the reversal of the effects of sense organs that signal forces on a leg when switching from forward to backward stepping in the stick insect (Akay et al., 2007, for review see Mantziaris et al., 2020). The mechanism underlying such a switch in sensory encoding could involve an interplay between the descending drive to the CPG and sensory feedback. In line with this possibility, brainstem stimulation changes how lamprey motoneurons respond to rhythmic movements imposed to the spinal cord (Hsu et al., 2013).

The regulation mechanism also explains the effect of limb feedback: the excitatory signal increases the frequency of the

limb oscillators. These in turn increase the frequency of the first segments, and thus the intersegmental phase lag. Such limb feedback has been proposed in a simulation study as a way of facilitating the transition from walking to trotting in the salamander (Harischandra et al., 2011).

The cells underlying proprioceptive axial feedback remain to be identified in salamanders (see section Motor control in salamanders). The limb sensory feedback introduced in simulation could be provided by cutaneous receptors during foot contact since mechano-sensitive Merkel cells are present on the skin of salamanders (Scott et al., 1981, Diamond et al., 1986), and/or stretch receptors of limb muscles that are sensitive to joint angle, since fibers behaving as muscle spindles have been identified in salamanders (Bone et al., 1976). In mammals, it is well-established that limb feedback plays a key role in establishing the locomotor patterns (e.g., Musienko et al., 2012; Akay et al., 2014; Takeoka et al., 2014; for review see Frigon, 2017).

## Muscles and Passive Biomechanical Properties

We found that higher muscle torques were required to emulate struggling and backward stepping (Table 5). Behavior-dependent changes in limb electromyographic activity have been reported in salamanders when comparing forward and backward terrestrial stepping. The electromyographic bursts increase during backward stepping in the extensor iliotibialis pars posterior (the homolog of the rectus femoris in mammals, which elevates the femur and extends the knee), mostly during the swing phase, whereas the bursts decrease in the other limb muscles (Ashley-Ross and Lauder, 1997). Future studies should determine whether an increase in electromyographic activity occurs in axial muscles during backward terrestrial stepping. A differential ratio of activation of epaxial vs. hypaxial muscles in the animal could also occur, as observed when comparing forward underwater stepping and forward terrestrial stepping in salamanders (Deban and Schilling, 2009). The same comparative electromyographic measurements should be done for struggling in salamanders. Caudorostral waves of axial activity are also used during struggling in Xenopus and during backward swimming in eels and lampreys. Lateral body undulations are much larger during struggling and backward swimming than during forward swimming in Xenopus (Kahn and Roberts, 1982), in eel (D'Aout and Aerts, 1999) and in lamprey despite a similar duty cycle of the electromyographic burst (Islam et al., 2006), suggesting that an increase in muscle strength occurs during caudorostral waves.

Passive tail segments reduced the drag during forward underwater stepping in simulation (Figure 7). In line with this, tail muscles show weak or no activation despite large tail undulations during forward underwater stepping in salamanders, suggesting that the body generates thrust by transmitting trunk movements passively to the tail (Cabelguen et al., 2014). At the low frequencies of underwater stepping, the passive biomechanical properties of the tail could be sufficient to propagate the body undulation, while higher frequencies might require a higher stiffness and thus active muscles (Blight, 1976, 1977) as observed in the salamander during swimming (Delvolvé et al., 1997). In salamanders, whether tail deactivation during forward underwater stepping is due to reduced activity of
e.g., some reticulospinal neurons remains to be determined. Biological observations and robotic experiments suggest that salamanders can use their tail as a “fifth limb” to provide thrust in slippery conditions (Karakasiliotis and Ijspeert, 2009). This suggests the existence of reticulospinal neurons that can both decrease and increase the activity of the tail independently from the trunk.

## Robotic Platform and Distributed Control

Some of the adaptations required to reproduce the simulation results on the robot can be explained by mechanical differences: The need for stronger muscle forces (no tapering) in the robot's tail might be due to the passive fin exerting more resistance than in the simulation. The different optimal limb-body phase bias for backward stepping could be due to the backlash in the leg gears.

Our initial implementation of the controller was centralized in the head module. This required retrieving the position and velocity of all joints and sending back the torque setpoints at each time step over the CAN bus. These are slow operations since the module has to forward the requests over a local $I^2C$ bus. The resulting control loop was too slow, making the muscle model unstable. The distributed controller solved this problem by keeping the communication of joint positions, velocities and torque setpoints local to each module. This solution shows interesting similarities with the vertebrate nervous system, which distributes the processing of sensory signals and the generation of locomotor patterns along the spinal cord, close to the target muscles.

Watanabe et al. (2009) have shown that a distributed controller with proprioceptive feedback can have interesting fault-tolerance properties, such as robustness to lesions in the communication pathways. It would be interesting to experiment with such lesions in our CPG model. The distributed controller would probably accommodate such experiments: the **Supplementary Movie 9** shows that the different sections of the CPG continue to function after the robot has been split in several parts. This is an interesting feature that few robots have. It is made possible by the distributed computation of the CPG and muscle model, the multi-master nature of the CAN bus and the nearest-neighbor couplings of the CPG model.

The distributed controller also introduces a difficulty in the form of coupling delays, which can be hard to predict when many modules share the same communication bus. As illustrated in **Figure 8A**, these delays can have a significant impact on the coordination between modules: the phase lags between rostral modules are markedly increased, while those between caudal modules are decreased. The asymmetry is probably related to the priority of messages on the CAN bus: the last modules have higher CAN identifiers so lower priorities when several modules attempt to talk at the same time. This means that ascending couplings will be on average more delayed than descending couplings, inducing larger lags in the rostral modules (**Figure 8B**). The problem was mostly solved by extrapolating in the receiver module the state of the oscillators at the origin of the couplings (see Results). Extrapolating these states is easily done in our CPG model, where the state variables are the phase and amplitude of the oscillators: during steady state locomotion, these variables, respectively, grow at an almost constant rate or

stay almost constant. It would be more difficult to extrapolate the state in a model without explicit phase variables. In conclusion, the CAN bus, being shared by all modules, limits the benefits of the distributed controller. A future revision of the robot should include direct communication between adjacent modules in addition to the shared bus, to fully realize the benefits of distributed control.

## CONCLUSION

Following the analogy proposed by Loeb (2001), our study suggests that the spinal cord is as a puppet on strings, and that a complex motor repertoire can be generated by pulling a limited set of “sensory” or “descending” strings, which in turn take advantage of a flexible spinal motor circuit.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## AUTHOR CONTRIBUTIONS

DR, JK, AI, and J-MC designed the study. JK and AC performed numerical modeling, simulations, and robotic experiments. JK, J-MC, AI, and DR analyzed data. JK, AI, and DR wrote the paper. All authors contributed to the article and approved the submitted version.

## FUNDING

The support from the European Community (LAMPETRA grant: FP7-ICT-2007-1-216100) and the Fondation pour la Recherche Médicale (FRM, DBC 20101021008) are acknowledged. DR receives support from the Canadian Institutes of Health Research (407083); the Fonds de la Recherche - Québec (FRQS Junior 1 awards 34920 and 36772); the Natural Sciences and Engineering Research Council of Canada (RGPIN-2017-05522 and RTI-2019-00628); the Centre de Recherche du Centre Hospitalier Universitaire de Sherbrooke (CHUS); the fonds Jean-Luc Mongrain de la fondation du CHUS, the Faculté de médecine et des sciences de la santé; the Centre d'excellence en Neurosciences de l'Université de Sherbrooke. AI receives support from Swiss National Science Foundation (Project CR23I2_140714).

## ACKNOWLEDGMENTS

The authors would like to thank Örjan Ekeberg for his valuable comments on a previous version of the manuscript, and Konstantinos Karakasiliotis for his help in preparing a movie of the robot. Part of this work was previously published as a thesis (Knüsel, 2013).

## SUPPLEMENTARY MATERIAL

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

Akay, T., Ludwar, B. C., Göritz, M. L., Schmitz, J., and Büschges, A. (2007). Segment specificity of load signal processing depends on walking direction in the stick insect leg muscle control system. *J. Neurosci.* 27, 3285-3294. doi: 10.1523/JNEUROSCI.5202-06.2007

Akay, T., Tourtelotte, W. G., Arber, S., and Jessell, T. M. (2014). Degradation of mouse locomotor pattern in the absence of proprioceptive sensory feedback. *Proc. Natl. Acad. Sci. U.S.A.* **111**, 16877-16882. doi: 10.1073/pnas.1419045111

Ampatzis, K., Song, J., Ausborn, J., and El Manira, A. (2014). Separate microcircuit modules of distinct v2a interneurons and motoneurons control the speed of locomotion. *Neuron* 83, 934-943. doi: 10.1016/j.neuron.2014.07.018

Ashley-Ross, M. A., and Lauder, G. V. (1997). Motor patterns and kinematics during backward walking in the pacific giant salamander: evidence for novel motor output. *J. Neurophysiol.* 78, 3047-3060. doi: 10.1152/jn.1997.78.6.3047

Ashley-Ross, M. A., Lundin, R., and Johnson, K. L. (2009). Kinematics of level terrestrial and underwater walking in the California newt, Taricha torosa. *J. Exp. Zool. A Ecol. Genet. Physiol.* 311, 240-257. doi: 10.1002/jez.522

Bejan, A., and Marden, J. H. (2006). Unifying constructal theory for scale effects in running, swimming and flying. *J. Exp. Biol.* 209, 238-248. doi: 10.1242/jeb.01974

Bem, T., Cabelguen, J., Ekeberg, Ö., and Grillner, S. (2003). From swimming to walking: a single basic network for two different behaviors. *Biol. Cybern.* 88, 79-90. doi: 10.1007/s00422-002-0340-3

Berg, E. M., Björnfors, E. R., Pallucchi, I., Picton, L. D., and El Manira, A. (2018). Principles governing locomotion in vertebrates: lessons from Zebrafish. *Front. Neural Circuits* 12:73. doi: 10.3389/fncir.2018.00073

Bicanski, A., Ryczko, D., Cabelguen, J.-M., and Ijspeert, A. J. (2013). From lamprey to salamander: an exploratory modeling study on the architecture of the spinal locomotor networks in the salamander. *Biol. Cybern.* 107, 565-587. doi: 10.1007/s00422-012-0538-y

Blight, A. (1976). Undulatory swimming with and without waves of contraction. *Nature* 264, 352-354. doi: 10.1038/264352a0

Blight, A. R. (1977). The muscular control of vertebrate swimming movements. *Biol. Rev.* 52, 181-218. doi: 10.1111/j.1469-185X.1977.tb01349.x

Böhm, U. L., Prendergast, A., Djenoune, L., Nunes Figueiredo, S., Gomez, J., Stokes, C., et al. (2016). CSF-contacting neurons regulate locomotion by relaying mechanical stimuli to spinal circuits. *Nat. Commun.* 7:10866. doi: 10.1038/ncomms10866

Bone, Q., Ridge, R. M., and Ryan, K. P. (1976). Stretch receptors in urodele limb muscles. *Cell Tissue Res.* 165, 249-266. doi: 10.1007/BF00226663

Branchereau, P., Rodriguez, J. J., Delvolvé, I., Abrous, D. N., Le Moal, M., and Cabelguen, J. M. (2000). Serotonergic systems in the spinal cord of the amphibian urodele pleurodeles waltl. *J. Comp. Neurol.* 419, 49-60. doi: 10.1002/(SICI)1096-9861(20000327)419:1&lt;49::AID-CNE3&gt;3.0.CO;2-#

Brocard, F., Ryczko, D., Fénelon, K., Hatem, R., Gonzales, D., Auclair, F., et al. (2010). The transformation of a unilateral locomotor command into a symmetrical bilateral activation in the brainstem. *J. Neurosci.* 30, 523-533. doi: 10.1523/JNEUROSCI.3433-09.2010

Buchanan, J. T. (2001). Contributions of identifiable neurons and neuron classes to lamprey vertebrate neurobiology. *Prog. Neurobiol.* 63, 441-466. doi: 10.1016/S0301-0082(00)00050-2

Buchanan, J. T., Grillner, S., Cullheim, S., and Risling, M. (1989). Identification of excitatory interneurons contributing to generation of locomotion in lamprey: structure, pharmacology, and function. *J. Neurophysiol.* 62, 59-69. doi: 10.1152/jn.1989.62.1.59

Büschges, A., Schmitz, J., and Bässler, U. (1995). Rhythmic patterns in the thoracic nerve cord of the stick insect induced by pilocarpine. *J. Exp. Biol.* 198, 435-456. doi: 10.1007/BF01185408

Cabelguen, J.-M., Bourcier-Lucas, C., and Dubuc, R. (2003). Bimodal locomotion elicited by electrical stimulation of the midbrain in the salamander *Notophthalmus viridescens*. *J. Neurosci.* 23, 2434-2439. doi: 10.1523/JNEUROSCI.23-06-02434.2003

Cabelguen, J.-M., Charrier, V., and Mathou, A. (2014). Modular functional organisation of the axial locomotor system in salamanders. *Zoology* 117, 57-63. doi: 10.1016/j.zool.2013.10.002

Caggiano, V., Leiras, R., Goñi-Erro, H., Masini, D., Bellardita, C., Bouvier, J., et al. (2018). Midbrain circuits that set locomotor speed and gait selection. *Nature* 553, 455-460. doi: 10.1038/nature25448

Cangiano, L., and Grillner, S. (2003). Fast and slow locomotor burst generation in the hemispheral cord of the lamprey. *J. Neurophysiol.* 89, 2931-2942. doi: 10.1152/jn.01100.2002

Cangiano, L., and Grillner, S. (2005). Mechanisms of rhythm generation in a spinal locomotor network deprived of crossed connections: the lamprey hemicord. *J. Neurosci.* 25, 923-935. doi: 10.1523/JNEUROSCI.2301-04.2005

Cangiano, L., Hill, R. H., and Grillner, S. (2012). The hemisegmental locomotor network revisited. *Neuroscience* 210, 33-37. doi: 10.1016/j.neuroscience.2012.03.007

Capelli, P., Pivetta, C., Soledad Esposito, M., and Arber, S. (2017). Locomotor speed control circuits in the caudal brainstem. *Nature* 551, 373-377. doi: 10.1038/nature24064

Charrier, V., and Cabelguen, J.-M. (2013). Fictive rhythmic motor patterns produced by the tail spinal cord in salamanders. *Neuroscience* 255, 191-202. doi: 10.1016/j.neuroscience.2013.10.020

Cheng, J., Stein, R. B., Jovanović, K., Yoshida, K., Bennett, D. J., and Han, Y. (1998). Identification, localization, and modulation of neural networks for walking in the mudpuppy (*Necturus maculatus*) spinal cord. *J. Neurosci.* 18, 4295-4304. doi: 10.1523/JNEUROSCI.18-11-04295.1998

Chevallier, S., Jan Ijspeert, A., Ryczko, D., Nagy, F., and Cabelguen, J.-M. (2008). Organisation of the spinal central pattern generators for locomotion in the salamander: biology and modelling. *Brain Res. Rev.* 57, 147-161. doi: 10.1016/j.brainresrev.2007.07.006

Cregg, J. M., Leiras, R., Montalant, A., Wanken, P., Wickersham, I. R., and Kiehn, O. (2020). Brainstem neurons that command mammalian locomotor asymmetries. *Nat. Neurosci.* 23, 730-740. doi: 10.1038/s41593-020-0633-7

Crespi, A., Karakasiliotis, K., Guignard, A., and Ijspeert, A. J. (2013). Salamandra robotica II: an amphibious robot to study salamander-like swimming and walking *gaits*. *IEEE Trans. Robot.* 29, 308-320. doi: 10.1109/TRO.2012.2234311

D'Aout, K., and Aerts, P. (1999). A kinematic comparison of forward and backward swimming in the eel anguilla anguilla. *J. Exp. Biol.* 202, 1511-1521.

Deban, S. M., and Schilling, N. (2009). Activity of trunk muscles during aquatic and terrestrial locomotion in *Ambystoma maculatum*. *J. Exp. Biol.* 212, 2949-2959. doi: 10.1242/jeb.032961

Delvolvé, I., Bem, T., and Cabelguen, J. M. (1997). Epaxial and limb muscle activity during swimming and terrestrial stepping in the adult newt, Pleurodeles waltl. *J. Neurophysiol.* 78, 638-650. doi: 10.1152/jn.1997.78.2.638

Delvolvé, I., Branchereau, P., Dubuc, R., and Cabelguen, J. M. (1999). Fictive rhythmic motor patterns induced by NMDA in an *in vitro* brain stem-spinal cord preparation from an adult urodele. *J. Neurophysiol.* 82, 1074-1077. doi: 10.1152/jn.1999.82.2.1074

Diamond, J., Holmes, M., and Nurse, C. A. (1986). Are merkel cell-neurite reciprocal synapses involved in the initiation of tactile responses in salamander skin? *J. Physiol.* 376, 101-120. doi: 10.1113/jphysiol.1986.sp016144

Ekeberg, Ö. (1993). A combined neuronal and mechanical model of fish swimming. *Biol. Cybern.* 69, 363-374. doi: 10.1007/BF01185408

Fagerstedt, P., Orlovsky, G. N., Deliagina, T. G., Grillner, S., and Ullén, F. (2001). Lateral turns in the lamprey. *II. activity of reticulospinal neurons during the generation of fictive turns*. *J. Neurophysiol.* 86, 2257-2265. doi: 10.1152/jn.2001.86.5.2257

Frigon, A. (2017). The neural control of interlimb coordination during mammalian locomotion. *J. Neurophysiol.* 117, 2224-2241. doi: 10.1152/jn.00978.2016

Frolich, L. M., and Biewener, A. A. (1992). Kinematic and electromyographic analysis of the functional role of the body axis during terrestrial and aquatic locomotion in the salamander Ambystoma tigrinum. *J. Exp. Biol.* 162, 107-130.

Gabriel, J. P., Ausborn, J., Ampatzis, K., Mahmood, R., Eklöf-Ljunggren, E., and El Manira, A. (2011). Principles governing recruitment of motoneurons during swimming in zebrafish. *Nat. Neurosci.* 14, 93-99. doi: 10.1038/nn.2704

Grillner, S. (1981). "Control of locomotion in bipeds, tetrapods, and fish," in *Handbook of Physiology. The Nervous System*, Motor Control, ed M. Brookhart (Bethesda, MD: American Physiological Society), 1179-1236. doi: 10.1002/cphy.cp010226
Grillner, S., and El Manira, A. (2020). Current principles of motor control, with special reference to vertebrate locomotion. *Physiol. Rev.* 100, 271–320. doi: 10.1152/physrev.00015.2019

Grillner, S., McClellan, A., and Sigvardt, K. (1982). Mechanosensitive neurons in the spinal cord of the lamprey. *Brain Res.* 235, 169–173. doi: 10.1016/0006-8993(82)90208-6

Grillner, S., Williams, T., and Lagerbäck, P. A. (1984). The edge cell, a possible intraspinal mechanoreceptor. *Science* 223, 500–503. doi: 10.1126/science.6691161

Hägglund, M., Dougherty, K. J., Borgius, L., Itohara, S., Iwasato, T., and Kiehn, O. (2013). Optogenetic dissection reveals multiple rhythmogenic modules underlying locomotion. *Proc. Natl. Acad. Sci. U.S.A.* 110, 11589–11594. doi: 10.1073/pnas.1304365110

Harischandra, N., Cabelguen, J.-M., and Ekeberg, Ö. (2010). A 3D musculo-mechanical model of the salamander for the study of different gaits and modes of locomotion. *Front. Neurorobot.* 4:112. doi: 10.3389/fnbot.2010.00112

Harischandra, N., Knuesel, J., Kozlov, A., Bicanski, A., Cabelguen, J.-M., Ijspeert, A., et al. (2011). Sensory feedback plays a significant role in generating walking gait and in gait transition in salamanders: a simulation study. *Front. Neurorobot.* 5:3. doi: 10.3389/fnbot.2011.00003

Harper, C. E., and Roberts, A. (1993). Spinal cord neuron classes in embryos of the smooth newt *Triturus vulgaris*: a horseradish peroxidase and immunocytochemical study. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* 340, 141–160. doi: 10.1098/rstb.1993.0053

Horvat, T., and Ijspeert, A. J. (2017). "Spine-limb coordination controller for a sprawling posture robot," in *8th International Symposium on Adaptive Motion of Animals and Machines (AMAM)* (Sapporo).

Horvat, T., Melo, K., and, and A. J., Ijspeert (2017). Spine controller for a sprawling posture robot. *IEEE Robot. Automat. Lett.* 2, 1195–1202. doi: 10.1109/LRA.2017.2664898

Hsu, L.-J., Zelenin, P. V., Grillner, S., Orlovsky, G. N., and Deliagina, T. G. (2013). Intraspinal stretch receptor neurons mediate different motor responses along the body in lamprey. *J. Comp. Neurol.* 521, 3847–3862. doi: 10.1002/cne.23382

Hubbard, J. M., Böhm, U. L., Prendergast, A., Tseng, P.-E. B., Newman, M., Stokes, C., et al. (2016). Intraspinal sensory neurons provide powerful inhibition to motor circuits ensuring postural control during locomotion. *Curr. Biol.* 26, 2841–2853. doi: 10.1016/j.cub.2016.08.026

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

Ijspeert, A. J., Crespi, A., and Cabelguen, J.-M. (2005). Simulation and robotics studies of salamander locomotion: applying neurobiological principles to the control of locomotion in robots. *Neuroinformatics* 3, 171–195. doi: 10.1385/NI:3:3:171

Ijspeert, A. J., Crespi, A., Ryczko, D., and Cabelguen, J.-M. (2007). From swimming to walking with a salamander robot driven by a spinal cord model. *Science* 315, 1416–1420. doi: 10.1126/science.1138353

Islam, S. S., Zelenin, P. V., Orlovsky, G. N., Grillner, S., and Deliagina, T. G. (2006). Pattern of motor coordination underlying backward swimming in the lamprey. *J. Neurophysiol.* 96, 451–460. doi: 10.1152/jn.0127.7.2005

Jalalvand, E., Robertson, B., Wallén, P., and Grillner, S. (2016). Ciliated neurons lining the central canal sense both fluid movement and pH through ASIC3. *Nat. Commun.* 7:10002. doi: 10.1038/ncomms10002

Josset, N., Roussel, M., Lemieux, M., Lafrance-Zoubga, D., Rastqar, A., and Bretzner, F. (2018). Distinct contributions of mesencephalic locomotor region nuclei to locomotor control in the freely behaving mouse. *Curr. Biol.* 28, 884–901.e3. doi: 10.1016/j.cub.2018.02.007

Kahn, J. A., and Roberts, A. (1982). The neuromuscular basis of rhythmic struggling movements in embryos of *Xenopus laevis*. *J. Exp. Biol.* 99, 197–205.

Karakasiliotis, K., and Ijspeert, A. J. (2009). "Analysis of the terrestrial locomotion of a salamander robot," in *2009 IEEE/RSJ International Conference on Intelligent Robots and Systems*. (St. Louis, MO), 5015–5020.

Karakasiliotis, K., Schilling, N., Cabelguen, J.-M., and Ijspeert, A. J. (2013). Where are we in understanding salamander locomotion: biological and robotic perspectives on kinematics. *Biol. Cybern.* 107, 529–544. doi: 10.1007/s00422-012-0540-4

Karakasiliotis, K., Thandiackal, R., Melo, K., Horvat, T., Mahabadi, N. K., Tsitkov, S., et al. (2016). From cineradiography to biorobots: an approach for designing robots to emulate and study animal locomotion. *J. R. Soc. Interface* 13:20151089. doi: 10.1098/rsif.2015.1089

Kiehn, O. (2016). Decoding the organization of spinal circuits that control locomotion. *Nat. Rev. Neurosci.* 17, 224–238. doi: 10.1038/nrn.2016.9

Kiemel, T., and Cohen, A. H. (2001). Bending the lamprey spinal cord causes a slowly-decaying increase in the frequency of fictive swimming. *Brain Res.* 900, 57–64. doi: 10.1016/S0006-8993(01)02244-2

Kinkhabwala, A., Riley, M., Koyama, M., Monen, J., Satou, C., Kimura, Y., et al. (2011). A structural and functional ground plan for neurons in the hindbrain of zebrafish. *Proc. Natl. Acad. Sci. U.S.A.* 108, 1164–1169. doi: 10.1073/pnas.1012185108

Knafo, S., Fidelin, K., Prendergast, A., Tseng, P.-E. B., Parrin, A., Dickey, C., et al. (2017). Mechanosensory neurons control the timing of spinal microcircuit selection during locomotion. *eLife* 6:e25260. doi: 10.7554/eLife.25260

Knüsel, J. (2013). *Modeling a Diversity of Salamander Motor Behaviors with Coupled Abstract Oscillators and a Robot*. Lausanne: EPFL.

Knüsel, J., Bicanski, A., Ryczko, D., Cabelguen, J.-M., and Ijspeert, A. J. (2013). A salamander's flexible spinal network for locomotion, modeled at two levels of abstraction. *Integr. Comp. Biol.* 53, 269–282. doi: 10.1093/icb/ict067

Kozlov, A., Huss, M., Lansner, A., Kotaleski, J. H., and Grillner, S. (2009). Simple cellular and network control principles govern complex patterns of motor behavior. *Proc. Natl. Acad. Sci. U.S.A.* 106, 20027–20032. doi: 10.1073/pnas.0906722106

Lavrov, I., and Cheng, J. (2004). Activation of NMDA receptors is required for the initiation and maintenance of walking-like activity in the mudpuppy (*Necturus Maculatus*). *Can. J. Physiol. Pharmacol.* 82, 637–644. doi: 10.1139/y04-044

Liu, Q., Yang, H., Zhang, J., and Wang, J. (2018). A new model of the spinal locomotor networks of a salamander and its properties. *Biol. Cybern.* 112, 369–385. doi: 10.1007/s00422-018-0759-9

Liu, Q., Zhang, Y., Wang, J., Yang, H., and Hong, L. (2020). Modeling of the neural mechanism underlying the terrestrial turning of the salamander. *Biol. Cybern.* doi: 10.1007/s00422-020-00821-1

Loeb, G. E. (2001). Learning from the spinal cord. *J. Physiol.* 533, 111–117. doi: 10.1111/j.1469-7793.2001.0111b.x

Lowry, C. A., Rose, J. D., and Moore, F. L. (1996). Corticotropin-releasing factor enhances locomotion and medullary neuronal firing in an amphibian. *Horm. Behav.* 30, 50–59. doi: 10.1006/hbeh.1996.0008

Mantziaris, C., Bockemühl, T., and Büschges, A. (2020). Central pattern generating networks in insect locomotion. *Dev. Neurobiol.* 80, 16–30. doi: 10.1002/dneu.22738

Matsushima, T., and Grillner, S. (1990). Intersegmental co-ordination of undulatory movements-a "trailing oscillator" hypothesis. *Neuroreport* 1, 97–100. doi: 10.1097/00001756-199010000-00003

Matsushima, T., and Grillner, S. (1992). Neural mechanisms of intersegmental coordination in lamprey: local excitability changes modify the phase coupling along the spinal cord. *J. Neurophysiol.* 67, 373–388. doi: 10.1152/jn.1992.67.2.373

McLean, D. L., Fan, J., Higashijima, S., Hale, M. E., and Fetcho, J. R. (2007). A topographic map of recruitment in spinal cord. *Nature* 446, 71–75. doi: 10.1038/nature05588

McLean, D. L., Masino, M. A., Koh, I. Y. Y., Lindquist, W. B., and Fetcho, J. R. (2008). Continuous shifts in the active set of spinal interneurons during changes in locomotor speed. *Nat. Neurosci.* 11, 1419–1429. doi: 10.1038/nn.2225

Mulloney, B., Harness, P. I., and Hall, W. M. (2006). Bursts of information: coordinating interneurons encode multiple parameters of a periodic motor pattern. *J. Neurophysiol.* 95, 850–861. doi: 10.1152/jn.00939.2005

Mulloney, B., and Smarandache-Wellmann, C. (2012). Neurobiology of the crustacean swimmeret system. *Prog. Neurobiol.* 96, 242–267. doi: 10.1016/j.pneurobio.2012.01.002

Musienko, P. E., Zelenin, P. V., Lyalka, V. F., Gerasimenko, Y. P., Orlovsky, G. N., and Deliagina, T. G. (2012). Spinal and supraspinal control of the direction of stepping during locomotion. *J. Neurosci.* 32, 17442–17453. doi: 10.1523/JNEUROSCI.3757-12.2012

Orts-Del'Immagine, A., Cantaut-Belarif, Y., Thouvenin, O., Roussel, J., Baskaran, A., Langui, D., et al. (2020). Sensory neurons contacting the cerebrospinal fluid
require the reissner fiber to detect spinal curvature *in vivo*. *Curr. Biol.* 30, 827-839. doi: 10.1016/j.cub.2019.12.071

Partridge, L. D. (1965). Modifications of neural output signals by muscles: a frequency response study. *J. Appl. Physiol.* 20, 150-156. doi: 10.1152/jappl.1965.20.1.150

Porez, M., Boyer, F., and Ijspeert, A. J. (2014). Improved Lighthill fish swimming model for bio-inspired robots: modeling, computational aspects and experimental comparisons. *Int. J. Robot. Res.* 33, 1322-1341. doi: 10.1177/0278364914525811

Ryczko, D., Auclair, F., Cabelguen, J.-M., and Dubuc, R. (2016a). The mesencephalic locomotor region sends a bilateral glutamatergic drive to hindbrain reticulospinal neurons in a tetrapod. *J. Comp. Neurol.* 524, 1361-1383. doi: 10.1002/cne.23911

Ryczko, D., Charrier, V., Ijspeert, A., and Cabelguen, J.-M. (2010). Segmental oscillators in axial motor circuits of the salamander: distribution and bursting mechanisms. *J. Neurophysiol.* 104, 2677-2692. doi: 10.1152/jn.0047 9.2010

Ryczko, D., and Dubuc, R. (2013). The multifunctional mesencephalic locomotor region. *Curr. Pharm. Des.* 19, 4448-4470. doi: 10.2174/1381612811319 240011

Ryczko, D., Knüsel, J., Crespi, A., Lamarque, S., Mathou, A., Ijspeert, A. J., et al. (2015). Flexibility of the axial central pattern generator network for locomotion in the salamander. *J. Neurophysiol.* 113, 1921-1940. doi: 10.1152/jn.0089 4.2014

Ryczko, D., Simon, A., and Ijspeert, A. J. (2020). Walking with salamanders: from molecules to biorobotics. *Trends Neurosci.* 43, 916-930. doi: 10.1016/j.tins.2020.08.006

Ryczko, D., Thandiackal, R., and Ijspeert, A. J. (2016b). "Interfacing a salamander brain with a salamander-like robot: Control of speed and direction with calcium signals from brainstem reticulospinal neurons," in *2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics* (Singapore: BioRob), 1140-1147.

Schilling, M., and Cruse, H. (2020). Decentralized control of insect walking: a simple neural network explains a wide range of behavioral and neurophysiological results. *PLoS Comput. Biol.* 16:e100780. doi: 10.1371/journal.pcbi.1007804

Schroeder, D. M., and Egar, M. W. (1990). Marginal neurons in the urodele spinal cord and the associated denticulate ligaments. *J. Comp. Neurol.* 301, 93-103. doi: 10.1002/cne.903010109

Scott, S. A., Cooper, E., and Diamond, J. (1981). Merkel cells as targets of the mechanosensory nerves in salamander skin. *Proc. R. Soc. Lond. B Biol. Sci.* 211, 455-470. doi: 10.1098/rspb.1981.0017

Smarandache, C., Hall, W. M., and Mulloney, B. (2009). Coordination of rhythmic motor activity by gradients of synaptic strength in a neural circuit that couples modular neural oscillators. *J. Neurosci.* 29, 9351-9360. doi: 10.1523/JNEUROSCI.1744-09.2009

Suzuki, S., Kano, T., Ijspeert, A. J., and Ishiguro, A. (2019a). "A salamander robot driven by cross-coupled sensory feedback control between legs and trunk," in *9th International Symposium on Adaptive Motion of Animals and Machines (AMAM)* (Lausanne).

Suzuki, S., Kano, T., Ijspeert, A. J., and Ishiguro, A. (2019b). Decentralized control with cross-coupled sensory feedback between body and limbs in sprawling locomotion. *Bioinspir. Biomim.* 14:066010. doi: 10.1088/1748-3190/ab3ef6

Takeoka, A., Vollenweider, I., Courtine, G., and Arber, S. (2014). Muscle spindle feedback directs locomotor recovery and circuit reorganization after spinal cord injury. *Cell* 159, 1626-1639. doi: 10.1016/j.cell.2014.11.019

Talpalar, A. E., Bouvier, J., Borgius, L., Fortin, G., Pierani, A., and Kiehn, O. (2013). Dual-mode operation of neuronal networks involved in left-right alternation. *Nature* 500, 85-88. doi: 10.1038/nature12286

Tytell, E. D., Hsu, C. Y., Williams, T. L., Cohen, A. H., and Fauci, L. J. (2010). Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming. *Proc. Natl. Acad. Sci. U.S.A.* 107, 19832-19837. doi: 10.1073/pnas.10115 64107

Watanabe, W., Sato, T., and Ishiguro, A. (2009). "A fully decentralized control of a serpentine robot based on the discrepancy between body, brain and environment," in *IEEE/RSJ International Conference on Intelligent Robots and Systems* (St. Louis, MO), 2421-2426.

Wyart, C., Del Bene, F., Warp, E., Scott, E. K., Trauner, D., Baier, H., et al. (2009). Optogenetic dissection of a behavioural module in the vertebrate spinal cord. *Nature* 461, 407-410. doi: 10.1038/nature08323

Yin, X., Jia, L., Wang, C., and Xie, G.-M. (2016). The effect of waist twisting on walking speed of an amphibious salamander like robot. *Acta Mech. Sin.* 32, 546-550. doi: 10.1007/s10409-015-0532-4

Zhong, B., Aydin, Y. O., Gong, C., Sartoretti, G., Wu, Y., Rieser, J., et al. (2018). Coordination of back bending and leg movements for quadrupedal locomotion. *Robot. Sci. Syst.* 14. doi: 10.15607/RSS.2018.XIV.020

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer CL declared a past co-authorship with one of the authors AI to the handling editor.

Copyright © 2020 Knüsel, Crespi, Cabelguen, Ijspeert and Ryczko. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.