# Robust undulatory locomotion through neuromechanical adjustments in a dissipative medium

Kenta Ishimoto, Clément Moreau, Johann Herault

► To cite this version:

Kenta Ishimoto, Clément Moreau, Johann Herault. Robust undulatory locomotion through neuromechanical adjustments in a dissipative medium. Journal of the Royal Society Interface, 2025, 22 (222), 10.1098/rsif.2024.0688. hal-05100551

HAL Id: hal-05100551

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

Submitted on 6 Jun 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
## INTERFACE

royalsocietypublishing.org/journal/rsif

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

Research

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

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

Check for
updates

**Cite this article:** Ishimoto K, Moreau C, Herault J. 2025 Robust undulatory locomotion through neuromechanical adjustments in a dissipative medium. *J. R. Soc. Interface* **22**: 20240688. https://doi.org/10.1098/rsif.2024.0688

Received: 2 May 2024
Accepted: 19 November 2024

**Subject Category:**
Life Sciences–Mathematics interface

**Subject Areas:**
biomechanics

**Keywords:**
animal locomotion, microswimmer, resistive force theory, coupled oscillators

**Authors for correspondence:**
Kenta Ishimoto
e-mail: ishimoto@kurims.kyoto-u.ac.jp
Clément Moreau
e-mail: clement.moreau@cnrs.fr
Johann Herault
e-mail: johann.herault@imt-atlantique.fr

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.7639319.

# Robust undulatory locomotion through neuromechanical adjustments in a dissipative medium

Kenta Ishimoto<fnref n="1" />, Clément Moreau<fnref n="2" /> and Johann Herault<fnref n="2" />

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

KI, 0000-0003-1900-7643; CM, 0000-0002-8557-1149; JH, 0000-0002-1536-7492

Dissipative environments are ubiquitous in nature, from microscopic swimmers in low-Reynolds-number fluids to macroscopic animals in frictional media. In this study, we consider a mathematical model of a slender elastic locomotor with an internal rhythmic neural pattern generator to examine various undulatory locomotion such as *Caenorhabditis elegans* swimming and crawling behaviours. By using local mechanical load as mechanosensory feedback, we have found that undulatory locomotion robustly emerges in different rheological media. This progressive behaviour is then characterized as a global attractor through dynamical systems analysis with a Poincaré section. Furthermore, by controlling the mechanosensation, we were able to design the dynamical systems to manoeuvre with progressive, reverse and turning motions as well as apparently random, complex behaviours, reminiscent of those experimentally observed in *C. elegans*. The mechanisms found in this study, together with our dynamical systems methodology, are useful for deciphering complex animal adaptive behaviours and designing robots capable of locomotion in a wide range of dissipative environments.

## 1. Introduction

Biological gait patterns vary in different media, such as swimming in water with fins, gliding and flying in air with wings and walking and crawling on the ground with legs. Among this morphological diversity, undulatory locomotion for swimming and crawling is used by various species, from flagellated cellular locomotion in invertebrates to swimming and crawling by vertebrates, such as fish and snakes [1–6].

It has long been a major challenge to find the common biological mechanisms that enable adaptation to different mechanical environments. In some vertebrates, such as lamprey, central pattern generators (CPGs) in the nervous system have been well known as the key mechanism underlying gait control and environmental adaptation. CPGs are made of networks of neural oscillators in the spinal cord, where autonomous neuronal circuits produce cyclic motor activity. This rhythmic pattern generates basic locomotion and can be modulated by sensory feedback to adapt to the environment [7–11].

Animal mechanosensation has been intensively studied in fish undulatory swimming at moderate and large Reynolds numbers, with a focus on proprioceptive feedback, which uses perceptions of local body deformation to sustain a desired gait pattern through stretch feedback [12–14]. Similarly, these aquatic animals can also perceive the dynamics of the fluid, giving rise to exteroceptive sensory feedback, through their lateral line or mechanical receptors [15–17]. The role of exteroceptive feedback in local sensorimotor loops is an ongoing research subject, but recent experiments with fish-like swimming robots in water [9] found that a robust undulatory locomotion can

© 2025 The Author(s). Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

<footnote n="1">¹Research Institute for Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="2">²Nantes Université, École Centrale Nantes, IMT Atlantique, CNRS, LS2N, UMR 6004, Nantes F-44000, France</footnote>
be achieved by local exteroceptive feedback control. These vertebrates live in a large-Reynolds-number regime, where the inertia of the body and fluid plays a crucial role. However, undulatory locomotion is used even when the inertia is negligible, and this regime is the focus of this study.

Animal locomotion in dissipative media includes swimming in low-Reynolds-number flow and crawling locomotion on a frictional surface or granular medium [18–20] (figure 1). In such a dissipation-dominated environment, drag is well described by the resistive force theory (RFT), which only considers local tangential and perpendicular drag coefficients. Here, we consider the resistive forces based on a viscous response of a Newtonian medium, as opposed to models based on Coulomb friction [23], such that the local force depends linearly on the velocity of the body section. This approach is known to be widely applicable across all scales in highly dissipative systems [19,20,24–27], with great accuracy for small deformations, that may decrease when non-local effects play an important role.

In the highly dissipative and inertialess regime, kinematic reversibility of the system necessitates non-reciprocal deformation of a self-propelled object to achieve net locomotion, which is well known as the scallop theorem [28,29] in low-Reynolds-number hydrodynamics. For instance, a complete in-phase synchronization of motor activation of ipsilateral muscles leads to a reciprocal deformation; hence, the animal should break the reciprocity of the internal states. In contrast, some fish can still swim in this case, a gait known as oscillatory swimming at large Reynolds numbers. Consequently, non-reciprocity in dissipative media imposes specific constraints on motor activity.

A strategy to overcome the reciprocity constraint would consist of coupling CPG dynamics with the fluid response to body motion through sensory feedback to produce elementary motor patterns for propulsion. Such a scenario could be exploited by microswimmers such as nematodes like *Caenorhabditis elegans*, although experimental evidence remains partial. Recent studies suggest the existence of a CPG-like network in *C. elegans*, reporting multiple sections of rhythmic cycles [30], and a relaxation-oscillator model supports these experimental findings [31].

The rhythmic neural activity of *C. elegans* has been used to develop integrated neuromechanical models for its undulatory locomotion [32,33]. More recently, based on RFT, fluid-structure interactions have been incorporated into models that contain both muscle activation and neural coupling with proprioceptive feedback. These models have successfully reproduced gait transitions of *C. elegans* between swimming and crawling and waveform modulation for swimming in media with different viscosities [34,35]. Johnson *et al.* [35] further theoretically derived a coarse-grained model from their integrated model by phase-oscillator reduction and found that the reduced model captures well the gait adaptation of swimming in *C. elegans*, motivating a simplified and minimal mathematical model for animal gait adaptation in dissipative media.

Hence, in this study, we develop a simple and widely applicable algorithm for neuromechanical adaption to different environments in a dissipation-dominated regime, by using local exteroceptive mechanical feedback from the environment. To do so, we extend the coupled-oscillator model with fluid-structure interactions used in [9], which reported robust self-organized locomotion through local hydrodynamic sensing in the inertial regime.

Our primary aim is to understand the mechanisms of neuromechanical adjustment and the complexity of animal behaviour in a dissipation-dominated regime through a simple and minimal mathematical model. Therefore, we introduce a high-level model of motor activation taking the form of CPGs that produce reciprocal motion in the absence of external stimuli. The underlying motivation is to introduce a fundamental motor control that does not favour any particular propulsive gait, due to motion reciprocity. By introducing sensory feedback in this specific CPG, we show that locomotion can be an emerging property of couplings between the body dynamics and the neuromotor system in dissipative media. Indeed, we report various motions, such as crawling and swimming at low and high viscosity, resulting from the adaptation to mechanical properties of the swimmer or drag of the dissipative medium.

To describe the diversity of the observed motions, we take inspiration from the behavioural-state space representation of *C. elegans*, recently introduced in an experimental study [36]. In this framework, its behaviour can be fully characterized by the state of its mechanical system and body shape. Our second aim is to analyse and design complex behaviours in terms of smooth trajectories produced by our neuromechanical model.

To reproduce behavioural variability, we hypothesize that this smooth trajectory is ruled by a non-autonomous dynamical system displaying different attractors. Here, the gait is represented by a periodic orbit or limit cycle. We show how couplings between body-environment dynamics and CPGs can generate transient dynamics between stable and unstable periodic orbits, which can reproduce the reversal of the motion direction. Motivated by the turning motion of *C. elegans*, we also introduce a phase-locking mechanism in the CPG to create a local fixed point in the state space corresponding to a posture producing the omega-turn manoeuvre [37]. Beyond the case of locomotion in *C. elegans*, our aim is to present a general theoretical framework to describe the emergence of complex motor patterns.

The content of the article is as follows: §2 introduces our mathematical model of a slender elastic object in a dissipative medium. The inner muscle activity is driven by the rhythmic pattern generated by CPGs, and the CPG phase is described by a coupled oscillator, which is locally modulated through mechanosensory feedback. We present typical emergent behaviours in §3 and further analyse the dynamical system with a focus on the trajectories in a Poincaré section in §4. In §5b, by controlling the phase dynamics, we evaluate a manoeuvring algorithm that reproduces omega turns of *C. elegans*. Concluding remarks are made in §6.

The Matlab code used for all numerical simulations conducted in this article is freely available [38].
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Examples of undulatory locomotion in a strongly dissipative medium. (a) Human sperm cells having whip-like motion without a neural network. Videomicroscopy from video 1 of the supplemental material of [18]. (b) Polychaeta worm, *Armandia brevis*, which can burrow in mud. Image adapted from [21] under CC BY-SA license (copyright 2024, Florida Museum of Natural History). (c) Swimming nematode captured and recorded by the authors at the north campus of Kyoto University. (d) Sandfish lizard that can swim in sand. Image adapted from [22] (copyright 2013, National Academy of Sciences).</figcaption></figure>

## 2. Model

### 2.1. Body–environment coupling

We model the animal's body as an inextensible elastic slender rod of length $L$. The centreline is denoted by its arclength $s \in [0, L]$, and its position by $\mathbf{x}(s, t)$. To calculate the force and torque balance for the local rod segment, contributions from elasticity, drag from the environment and inner activation are incorporated.

We assume planar locomotion (we set this as our $xy$ plane) and no twisting motion. In biological locomotors, muscular activation usually occurs through the actuation of several elements along the body. For instance, C. elegans has approximately six muscular modules on its body [35]. This discrete activation structure within a continuous elastic rod is conveniently described by the coarse-grained representation of the elastic slender rod, known as the $N$-link model, where the object is represented by $N$ equal-length links connected at $N-1$ hinges [39,40] (figure 2). Here, we take $N = 10$, following a previous study [9]. We write the position of the end of the first link as $(X, Y)$ and the angle from the $x$-axis as $\theta$. As shown in figure 2, these represent the position and orientation of the body and define the body-fixed frame $\{\mathbf{e}_{x0}, \mathbf{e}_{y0}, \mathbf{e}_{z0}\}$, which is distinct from the laboratory frame $\{\mathbf{e}_{x}, \mathbf{e}_{y}, \mathbf{e}_{z}\}$. The object shape is described by $N-1$ relative angles between the neighbouring links and denoted by $\alpha_i$ ($i = 1, 2, ..., N-1$).

For the body-environment coupling, we use the RFT, in which the local force per unit length at $\mathbf{x}$ is assumed to be proportional to the local tangent and normal velocities, $\mathbf{u}_{\parallel}$ and $\mathbf{u}_{\perp}$, as $\mathbf{f}^{\text{env}} = -c_{\parallel}\mathbf{u}_{\parallel} - c_{\perp}\mathbf{u}_{\perp}$, where the coefficients, $c_{\parallel}$ and $c_{\perp}$, are the drag coefficient in each direction (figure 2 inset). Since the drag force by the environment is linearly dependent on the local velocity, to represent the local velocity, we need linear and angular velocities as well as deformation velocity, for which we introduce a generalized velocity of the moving and deforming body as $\dot{\mathbf{z}} = (\dot{X}_0, \dot{Y}_0, \dot{\theta}, \dot{\alpha}_1, ..., \dot{\alpha}_{N-1})^{\text{T}}$, where the dot symbol indicates the time derivative and $(\dot{X}_0, \dot{Y}_0)$ is the velocity vector in the body-fixed frame.

For inertialess motion, the dynamical system takes the form of a linear system between the generalized velocities and the internal torques, thanks to the generalized resistance matrix $\mathbf{A}(\boldsymbol{\alpha})$ [41]. This linear relationship is based on the proportionality between drag and generalized velocity, with coefficients depending on the instantaneous shape $\boldsymbol{\alpha}$. The matrix $\mathbf{A}(\boldsymbol{\alpha})$ is obtained by imposing balance of forces and torques, which includes drag and internal torques, on the whole body as well as on each junction, thus, giving $N+2$ linear equations coupling the generalized velocities to the internal torques. Detailed expressions are provided in [39]. We then arrive at a matrix form of the body–environment coupling as [42,43]:

$$ \dot{\mathbf{A}}\dot{\mathbf{z}} = \mathbf{t} \qquad (2.1) $$

The vector $\mathbf{t}$ on the right-hand side is the generalized force and torque vector, and for free locomotion only the bottom $N-1$ rows may have non-zero values as $\mathbf{t} = (0, 0, 0, \tau_1, ..., \tau_{N-1})^{\text{T}}$. Here, the torque at the $i$th hinge contains elastic bending and internal actuation, which we denote $\tau_i = \tau_i^{\text{ela}} + \tau_i^{\text{int}}$.

### 2.2. Coupled oscillators and mechanosensation

We assume that the inner activation is controlled by a network of interconnected neural oscillators forming the CPGs, which provide rhythmic motions to the junctions. To model the cyclic activation of the junctions, we assume that the $i$th oscillator is associated with the corresponding hinge, and generates a periodic activation signal given by a phase denoted as $\phi_i(t)$ ($i = 1, 2, ..., N-1$) (figure 2).

To model the elasticity of the body, we use the Euler–Bernoulli constitutive relation, which assumes a linear relation between the bending moment and local curvature. With a bending modulus $\kappa$ for the discretized rod, we then have $\tau_i^{\text{ela}} = -\kappa\alpha_i$. For torque actuation, we simply assume a periodic sinusoidal mode [9], given by $\tau_i^{\text{int}} = \tau \cos \phi_i(t)$, where $\tau$ is the torque actuation strength at each hinge. The torque function is then summarized as $\tau_i = -\kappa\alpha_i + \tau \cos \phi_i(t)$.

Following previous studies of CPGs for undulatory locomotion [44,45] to describe this system, we describe the oscillator dynamics by a limit cycle close to a Hopf bifurcation. In this case, phase reduction of the dynamics leads to [46]

$$ \dot{\phi}_i = \omega_0 + \sum_j C_{ij} \sin(\phi_j - \phi_i) + Z(\phi_i) S_i(t), \qquad (2.2) $$
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Schematic of the model and the two frames of reference. A slender elastic rod is represented by N links connected by elastic hinges. Each hinge generates internal torque $\tau_i(i = 1,2, ..., N-1)$, as a function of the phase $\phi_i(t)$, which models the local CPGs. The phase dynamics are described by a coupled-oscillator equation with modulation by local signal $S_i(t) = M_i^{(i)} + M_{i+1}^{(i)}$, with $M_i^{(j)}$ being the drag on the jth link around the jth hinge (inset). The laboratory-fixed frame $\{e_x, e_y, e_z\}$ and body-fixed frame $\{e_{x_0}, e_{y_0}, e_{z_0}\}$ are introduced, and the motions are assumed to be restricted in the x-y plane. The motion of the object is represented by the end of the first link (X, Y) and its angle to the laboratory frame $\Theta$. The shape of the body is described by $N-1$ relative angles $\alpha_i$.</figcaption></figure>

where we consider a symmetric, nearest-neighbour phase coupling with $C_{i,i \pm 1} = C$ and zero otherwise. The phase oscillators are assumed to have a common intrinsic frequency $\omega_0$ and coupling with their neighbours with a common strength that is modulated by sensory feedback from the surrounding mechanical environment. Close to a Hopf bifurcation, the phase sensitivity function takes the form of $Z(\phi_i) = \sigma \cos \phi_i(t)$ with a sensitivity strength $\sigma$ [9].

The phase sensitivity function $Z(\phi_i)$ characterizes the response of a limit cycle weakly perturbed by the sensory feedback $S_i(t)$. The function $Z(\phi_i)$ modulates the input signal $S_i(t)$ such that their correlation modifies the phase dynamics, resulting in-phase acceleration or deceleration. Synchronization and phase locking between these two signals is an emergent phenomenon that depends on both the choice of the sensitivity function and the nature of the sensory feedback.

As for the sensory feedback, we follow previous literature on biological exteroception [9,47,48]. In their model for rhythmic contraction of true slime mould, Kobayashi et al. [47,48] considered local peristaltic pumping generated by phase oscillators and used the local hydrodynamic pressure as exteroceptive feedback. In the study of Tandiackal et al. [9], in which exteroceptive feedback in fish-like robot is studied, the authors consider local hydrodynamic forces at the actuation point on both sides of the elongated body and use the difference between the forces on each side as the exteroceptive signal. In our model, taking the slender limit of this difference, we employ the local torque load at each hinge as the exteroceptive feedback. Let $M_i^{(j)}$ be torque (per unit length) from the environmental drag on the ith link around jth hinge, or equivalently,

$$M_i^{(j)} = e_z \cdot \int_{S_{i-1}}^{S_i} (\mathbf{x}(s) - \mathbf{x}(s_j)) \times \mathbf{f}^{\text{env}}(s) \frac{ds}{\Delta s'} \quad (2.3)$$

where $\Delta s = L/N$ and $S_i = i\Delta s$ corresponds to the position of the ith hinge (inset of figure 2). We then set a fore-aft symmetric feedback $S_i(t) = M_i^{(i)} + M_{i+1}^{(i)}$, which is the local torque load on the ith hinge.

Other choices to model proprioceptive sensory feedback include using the local curvature as signal input, which, however, must be asymmetric to break the fore-aft symmetry and generate undulatory locomotion. Some mathematical models consider asymmetric and non-local signal inputs with weight on the posterior neighbour (e.g. [35]), enforcing a directed motion. In contrast, our exteroceptive feedback is a local and symmetric function, which fits our purpose to simplify the mechanism for adaptive locomotion.

## 2.3. Non-dimensionalization and model parameters

The mechanical torque on the body is calculated through the RFT, which is applied for various undulatory motions in dissipation-dominated environments, such as low-Reynolds-number flows, gel-like structures and mud and sand. We then introduce the anisotropic drag ratio as $\gamma = c_{\parallel}/c_{\perp}$, which is known to affect locomotion. Following the data obtained for C. elegans, we use $\gamma = 1/2$ for swimming and $\gamma = 1/70$ for crawling [25,26]. When the model swims in a low-Reynolds-number Newtonian fluid, the coefficients are theoretically estimated [49] as $c_{\perp} = 4\pi\mu/[\log(2L/d)]$ in a slender asymptote, where $\mu$ is the fluid viscosity and $d$ is the cross-sectional radius of a slender body.

For our numerical and theoretical analyses, we non-dimensionalize the above system by setting the scales $L = 1$ for length, $T = 2\pi/\omega_0 = 1$ for time, and $\kappa = 1$ for force. An important non-dimensional number of the system is the so-called 'sperm number',
which represents the effective flexibility of an oscillatory elastic rod in a highly damped system, defined as $Sp = L(\omega_0 c_1 / \kappa_b)^{1/4}$ [39], where $\kappa_b = \kappa / (L/N)$ is the bending modulus. For flagellated microswimmer, a typical value of Sp for human sperm flagella in water is $Sp \approx 4$ [50]. The estimation of Sp for *C. elegans* locomotion; however, is limited due to estimated values of the bending modulus, which could vary by up to two orders of magnitude depending on experiments [35].

Sznitzman et al. [51] estimated the bending modulus as $\kappa_b \approx 4.2 \times 10^{-16} \text{N} \cdot \text{m}^2$ and observed the beating frequency $1/T \approx 2.4$ Hz in watery medium ($\mu \approx 10^{-3} \text{Pa} \cdot \text{s}$). With these parameters and $L = 10^{-3} \text{m}$ for the body length, we may estimate the value of sperm number as $Sp \approx 2.4$. Fang-Yen et al. [52], however, estimated the bending modulus as $\kappa_b \approx 9.5 \times 10^{-14} \text{N} \cdot \text{m}^2$ observe the beating frequency $1/T \approx 1.8$Hz. With these parameters, the sperm number is estimated as $Sp \approx 0.7$. Fang-Yen et al. [52] also examined *C. elegans* swimming in viscous medium. From their observations, one may estimate the value of sperm numbers as $Sp \approx 5.2$ for the medium $10^3$-fold more viscous than water with $1/T = 0.5$Hz, and as $Sp \approx 12.5$ for the $10^4$-fold viscous medium with $1/T = 0.15$Hz. Shen et al. [53] evaluated the drag coefficient of *C. elegans* crawling on agar gels, suggesting an increase of drag coefficient by circa 400 times than that of swimming in water. We may then estimate the sperm numbers for crawling ranging $Sp = 4.1 - 16.0$ by considering the decrease of the beat frequency in crawling ($1/T \approx 0.5$Hz). In this article, to cover a biologically relevant range, we consider a range of $Sp = 1 - 12$. Other parameters are drag coefficient ratio $\gamma$, non-dimensionalized actuation strength $\hat{\tau}$, coupling strength $\hat{C}$ and sensitivity strength $\hat{\sigma}$ using the same physical units of $L, T$ and $\kappa$. Hereafter, we omit the hat symbol for the non-dimensionalized strengths.

# 3. Robust undulatory locomotion

## 3.1. Breakdown of motion reciprocity

Before proceeding to the simulation results, we first consider how motion non-reciprocity can be produced by breaking the symmetry with respect to the locomotor's transverse plane in the oscillator chain from the fore-aft symmetric sensory feedback function. The schematic in figure 3 explains the symmetry-breaking mechanism. In figure 3a, the oscillators are fully synchronized with $\cos \phi_i > 0$ such that the anticlockwise torque generates a concave shape, where the end links rotate more quickly. The local torque load is therefore anticlockwise at the left end ($S_i > 0$), while the torque load is clockwise at the right end ($S_i < 0$), yielding a negative gradient of resistive torque along the rod. This gradient is then projected along the oscillator through the correlation $\cos(\phi)S_i$ in the phase equation (2.2).

When the synchronized phases are satisfied $\cos \phi_i < 0$ (figure 3b), the direction of the local torque actuation reverses and the shape becomes convex. As in figure 3b, the local torque load has the opposite sign in the left and right ends but with the opposite sign of $S_i$ for the case of figure 3a, that is, a positive gradient of $S_i$. Nonetheless, the modulation to the phase dynamics, a product of $S_i$ and $\cos \phi_i$, remains positive at the left end and negative at the right end in both phases. Thus, when the sensitivity strength, $\sigma$, is positive (negative), the phases are accelerated on the left (right) side of the object, and decelerated on the right (left) side of the object, leading to symmetry breaking at the origin of the travelling wave of motor activation.

## 3.2. Robust undulatory locomotion

We now present simulation results of our model, equations (2.1) and (2.2). We have numerically calculated by an explicit time integration with the variable order method, implemented by a MATLAB built-in solver (ode15s).

We evaluated the simulation results in a different medium by changing the sperm number $Sp$ and torque actuation strength $\tau$ for swimming dynamics with $\gamma = 1/2$. In a wide parameter regime, we found that the object eventually exhibits stable undulatory locomotion by breaking the fore-aft symmetry (figure 3). The direction of motion is governed by the sensitivity strength sign, and we found forward motion (moving towards the first link) when $\sigma > 0$ and backward motion (towards the last link) when $\sigma < 0$.

Obtained shape morphology and motion parameters of emergent undulatory locomotion are summarized in figure 4. Other model parameters are fixed as $C = 1$, $\sigma = 8$.

We present mean maximum angle between links, $A$, mean velocity, $V_p$, motion frequency, $T_p$, wavelength, $\lambda$ and head-to-tail phase difference $\Delta\phi$. Each parameter is schematically shown in the inset of figure 4d.

As shown in figure 4a, the values of $A$, which measures the maximum curvature and amplitude of the undulation, increases as $\tau$ increases, while these values decrease at a higher $Sp$ because of a larger viscous drag force. Some data at a small $Sp$ are not obtained due to non-physical large amplitude with collision between links.

The mean velocity decreases at a higher $Sp$ and drops to almost zero at a value of $Sp$ ranging between 6 and 10 depending on $\tau$ (figure 4b). These critical values of $Sp$ are clearly shown by a sudden drop of wavelength shown by arrows in figure 4c. Above these critical values of $Sp$, the coherent waves cease and locomotive behaviours are no longer observed. These changes of regime are reflected in saturation of $T_p$ around $T_p \approx 1$ in figure 4c and sudden change of the phase difference, $|\Delta\phi|$, in figure 4e.

For a fixed internal actuation $\tau$, we have found there is a wide range of $Sp$ with stable undulatory swimming but with a maximum $Sp$. To swim in a high $Sp$ region, one needs a larger actuation. Interestingly, this tendency of larger actuation in higher viscosity is compatible with *C. elegans* swimming [52] and human sperm swimming [54,55]. Moreover, the wavelength of undulatory swimming at a higher viscosity is known to decrease in these species, which also qualitatively agrees with the current results.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Schematic mechanism of symmetry breaking through mechanosensory feedback. Fully synchronized situations are illustrated with (a) cos φᵢ > 0 and (b) cos φᵢ < 0. The direction of the local torque actuation is opposite, and the shape becomes concave in (a) and convex in (b). Despite the uniform distribution of the torque, the local torque load is no more uniform along the object, yielding fore-aft asymmetric phase modulation due to mechanoreception. In both phases, the product of Sᵢ and cos φᵢ is always positive at the left end and negative at the right end. Thus, when the sensitivity strength, σ, is positive, the phases are accelerated at the left and decelerated at the right, leading to a travelling undulatory motion.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Shape morphology and motion parameters of emergent undulatory locomotion in different Sp with 1 ≤ Sp ≤ 12 and actuation strength τ ∈ {4, 8, 12, 16}. (a) Mean maximum angle, A, (b) mean velocity, *V*<sub>p</sub>, (c) motion frequency, *T*<sub>p</sub>, (d) wavelength, λ and (e) head-to-tail phase difference Δφ. Schematics of the parameters are shown in the inset of (d). With a small Sp and larger actuation, data are not shown due to non-physical large amplitude with body's self-intersection. The arrows in (d) show critical values of Sp that exhibit sudden drops of λ.</figcaption></figure>

For constant wavelength, we should observe |Δφ| = 2π/λ, which is almost true for small Sp with coherent undulatory waves. However, this does not hold when the coherency disappears for larger Sp, since the net phase lag no longer reflects the wavelength of body deformation.
The shape morphology and motion parameters for crawling motions with $\gamma = 1/70$ are qualitatively the same as the swimming motion of figure 4 with a motion decoherency at a large Sp accompanied by decrease of wavelength (see electronic supplementary material, figure S1). The dependence of other parameters such as C and $\sigma$, and $\gamma$ are also shown in electronic supplementary material, figures S2–S4, showing robust emergence of undulatory motion in a wide parameter regime. In particular, decreasing $\gamma$ in general leads to a larger undulatory motion wavelength in crawling (electronic supplementary material, figure S4).

In figure 5a, we show typical swimming dynamics ($\gamma = 1/2$) in a fluid with a lower (Sp = 4) and higher (Sp = 8) sperm numbers (see also electronic supplementary material, Videos S1-S2). The top panel shows snapshots of the shape of the swimmer in different colours. The end of the first link in each snapshot is denoted by a black dot for illustrative purposes. The model animal, initially aligned in a straight configuration, starts to create a travelling wave and eventually exhibits stable periodic swimming. The middle and bottom panels present the kymographs of the angle $\alpha_i(t)$ and the phase $\phi_i(t)$. The other parameters were set as (a) $\tau = 8$, $\sigma = 4$ and $C = 1$ and (b) $\tau = 16$, $\sigma = 12$.

We then focus on crawling behaviour by setting the drag anisotropy ratio to $\gamma = 1/70$ [26]. Example crawling motions are shown in figure 6 and electronic supplementary material, videos S3–S4. Of particular note, the locomotion trajectory begins to follow the body shape, which is seen in C. elegans [56]. In figure 6, we plot superimposed snapshots in the physical x-y space to show the shape and the trajectory of the end of the first link, which is again denoted by a black dot in each snapshot. The parameter set used in the figure was Sp = 8, $\tau = 16$, $\sigma = 12$, $C = 1$ in (a) and Sp = 6, $\tau = 8$, $\sigma = 12$, $C = 0.1$ in (b).

## 4. Stability of periodic orbits

In the previous section, stable, robust periodic locomotion was observed in a wide range of parameters. We now seek mathematical characterization of the existence and stability of these emerging periodic orbits, by adopting a dynamical system viewpoint.

We introduce the dynamics state for a generic N-link model as

$$ \mathbf{x} = (\alpha_1, \alpha_2, \dots, \alpha_{N-1}, \phi_1, \phi_2, \dots, \phi_{N-1})^T \in U, \quad (4.1) $$

where $U = \mathbb{R}^{N-1} \times \mathbb{T}^{N-1}$ is the state space. We choose our Poincaré section at $\phi_1 = 0$ and denote the section as $P = \mathbb{R}^{N-1} \times \mathbb{T}^{N-2}$, which leads to the associated Poincaré map $F: P \to P$ (figure 7b). Here, we consider a periodic orbit in the state space including its configuration and CPG phases, while the CPG phase is a coarse-grained representation of the oscillatory neural dynamics at each junction.

A fixed point $\mathbf{x}^*$ on the Poincaré section corresponds to a periodic orbit in the state space, or in other words, the periodic animal morphology, including undulatory locomotion. To detect periodic orbits in $U$, which is regarded as a fixed point in $P$ such that $F(\mathbf{x}^*) = \mathbf{x}^*$, we use the standard Newton-Raphson method through $\mathbf{x}_{m+1} = \mathbf{x}_m - [\mathbf{J}(\mathbf{x}_m)]^{-1}[F(\mathbf{x}_m) - \mathbf{x}_m]$ for an integer $m$, where $\mathbf{J}(\mathbf{x}_m)$ represents the Jacobian matrix. The eigenvalues of the Jacobian matrix approximate the Floquet exponents when the algorithm converges by iteration. For methodological details, one may refer to a review paper [57] and applications in sperm swimming dynamics [58].

### 4.1. Purcell's three-link model

We first consider the case of $N = 3$, where the three links are connected by two hinges (figure 7a). This is Purcell's three-link model [28], which has been widely used as a simple and canonical mathematical model of microscale swimming [42,59]. The state space $U$ is then four-dimensional and the Poincaré section $P$ is three-dimensional.

Before proceeding to the numerical results of the periodic orbits, let us first analyse the condition under which the periodic orbits exist. When $N = 3$, the phase dynamics are written as

$$ \left\{ \begin{aligned} \dot{\phi}_1 &= \omega_0 + C \sin(\phi_2 - \phi_1) + \sigma \cos(\phi_1)S_1(t) \\ \dot{\phi}_2 &= \omega_0 + C \sin(\phi_1 - \phi_2) + \sigma \cos(\phi_2)S_2(t)' \end{aligned} \right. \quad (4.2) $$

and by introducing the phase difference $\Delta\phi = \phi_2 - \phi_1$, we obtain the phase equation:

$$ \Delta\dot{\phi} = -2C \sin(\Delta\phi) + \sigma[\cos(\phi_2)S_2 - \cos(\phi_1)S_1]. \quad (4.3) $$

In the absence of mechanosensory feedback, namely, $\sigma = 0$, a complete in-phase synchronization, $\Delta\phi = 2n\pi$ ($n = 0, \pm 1, \pm 2, \dots$), is a stable fixed point, while a complete anti-phase synchronization, $\Delta\phi = (2n - 1)\pi$, is an unstable fixed point.

When the mechanosensory feedback is activated ($\sigma \neq 0$), two synchronized states can still be observed, which correspond to the in-phase and anti-phase synchronized states when $\sigma = 0$, and they inherit their respective stability properties. In the small-amplitude case ($\tau \ll 1$), we may analytically prove that these two synchronized states always exist, regardless of the values of Sp, $\gamma$, $\tau$, $\sigma$ and C. The detailed derivations are provided in the electronic supplementary materials.

In the finite-amplitude case, however, the synchronized states can cease due to a saddle-node bifurcation. Indeed, there may be a critical coupling constant $C^*$, below which the synchronized states disappear. Figure 8a plots the values of the critical coupling constant for different sizes of $\tau$ with the other parameters being the same as in figure 5a. At a large amplitude with $\tau \gtrsim 5.8$, phase slip occurs in the parameter region without synchronized states.
<figure><img src="image_10.png" /></figure>

**Figure 5.** Typical emergent motions of (a) swimming in a lower viscosity medium (Sp=4) and (b) swimming in a higher viscosity medium (Sp=8). The model animal was initially aligned in a straight configuration with the same phase. The top panels show the time series of shape gaits. The shape of the locomotor is plotted in different colours to present the time evolution and the black dot denotes the end of the first link in each snapshot. The time series of the angle $\alpha_i$ and the phase $\phi_i$ are shown as kymographs in the middle and bottom panels. The model parameters are (a) Sp = 4, $\tau$ = 8, $\sigma$ = 4, C = 1 and (b) Sp = 8, $\tau$ = 16, $\sigma$ = 12, C = 1. The resulting motion parameters of undulatory locomotion are (a) $A$ = 0.64, $\lambda$ = 1.67, $V_p$ = 0.109, $T_p$ = 0.952, $\Delta\phi/(2\pi)$ = -0.63, (b) $A$ = 0.95, $\lambda$ = 0.82, $V_p$ = 0.065, $T_p$ = 0.907, $\Delta\phi/(2\pi)$ = -1.31. See also electronic supplementary material, Videos S1-S2.

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

**Figure 6.** Typical emergent motions of crawling on an agar gel. The model animal was initially aligned in a straight configuration with the same phase. Snapshots of the shape in the xy plane are superimposed on the trajectory of the end of the first link. The shape of the locomotor is plotted in different colours to present the time evolution and the black dot denotes the end of the first link in each snapshot. The model parameters are (a) Sp = 8, $\tau$ = 16, $\sigma$ = 12, C = 1 and (b) Sp = 6, $\tau$ = 8, $\sigma$ = 12, C = 0.1. The resulting motion parameters of undulatory locomotion are (a) $A$ = 1.16, $\lambda$ = 1.25, $V_p$ = 0.288, $T_p$ = 0.892, $\Delta\phi/(2\pi)$ = -0.78, (b) $A$ = 0.76, $\lambda$ = 0.97, $V_p$ = 0.501, $T_p$ = 0.915, $\Delta\phi/(2\pi)$ = -1.08. See also electronic supplementary material, Videos S3-S4.

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

**Figure 7.** (a) Schematic of Purcell's three-link model. (b) Schematic of Poincaré section *P* in the state space *U* and the associated Poincare map *F*: *P* → *P*.

Nonetheless, this weak coupling region is very limited. In a large parameter region, which specifically contains the typical parameters for biologically relevant motion, there is one stable periodic orbit (SPO) and one unstable periodic orbit (UPO). We established through numerical stability analysis that the UPO is a saddle point. Figure 8b shows a schematic of the structure of the dynamical systems in *U* and on *P*, where the SPO (the red point) and UPO (the blue point) are connected by heteroclinic orbits. By setting the parameter values as those of the low-viscosity swimmer in figure 5a, we numerically detect the SPO and UPO. Figure 8c plots them in the $\alpha_1 - \alpha_2$ shape space, with red and blue orbits representing the SPO and UPO, respectively. When the initial shape is located near the UPO, the dynamics evolve along a heteroclinic orbit, as shown in figure 8d. As time progresses (from blue to red), the trajectory departs from the UPO and finally approaches the SPO with the cyclic direction changing from anticlockwise to clockwise, to anticlockwise in the $\alpha_1 - \alpha_2$ shape space. Time evolution of shapes and trajectories in real space for figure 8d are found in electronic supplementary material, video S5.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** (a) Plots of critical value $C^*$ for different $\tau$ values for Purcell's swimmer in a low-viscosity medium. (b) Schematic of the dynamical systems on Poincaré section *P*. There is one stable fixed point (red) and one saddle point (blue). (c) Periodic orbits of the three-link swimmers in the $\alpha_1 - \alpha_2$ shape space. The large orbit (red) is stable, and the thinner orbit (blue) is unstable. Arrows indicate the direction of the trajectories. (d) Time evolution in the $\alpha_1 - \alpha_2$ shape space, with the initial shape being set close to the unstable saddle. Different colours illustrate the time evolution. See also electronic supplementary material, video S5.</figcaption></figure>

## 4.2. Stable and unstable periodic orbits in N-link model: an example

We now proceed to the analysis of the generic N-link system with $N > 3$. In this subsection, to illustrate the dynamical systems near a saddle-node bifurcation, we set up $Sp = 6$ and $\tau = 8$ with a small $C/\sigma$ value as $C = 0.1$ and $\sigma = 12$ with $N = 10$. The emergent stable undulatory motions are shown in figure 6b and electronic supplementary material, video S4.

With the same numerical analysis of the Poincaré section, we detected a global SPO, by confirming that all the Floquet exponents have a negative real part. A snapshot of the shape gait is shown in the top left of figure 9a. In the same parameter regions, we also found a UPO. As is shown in the top right of figure 9a, the snapshot for this UPO shows a longer wavelength ($\lambda = 1.27$) than that for the SPO waveform ($\lambda = 0.97$). Other motion parameters are obtained as $A = 0.61$, $V_p = 0.667$, $T_p = 0.938$, $\Delta\phi/(2\pi) = -0.72$. We confirm that one of the Floquet exponents is a positive real number, indicating that the UPO is a saddle, like in the Purcell three-link model.

Hence, if the model configuration is initially close to the UPO, the shape and phase dynamics evolve along the heteroclinic orbits that connect the SPO and UPO, as shown in figure 9. The waveform shifts from the unstable longer wave to the stable shorter wave. See also electronic supplementary material, video S6.

In the parameter set of figure 9, the critical coupling constant, below which synchronization ceases, is observed to be very small, at $C^*/\sigma \approx 0.01$. Furthermore, we found that large-scale stable locomotion emerges even for a larger $N$ ($\approx 100$), suggesting that only one SPO exists even for a larger $N$.

In this section, we have seen a global attractor in the behavioural-state space, together with transient behaviour. In the next section, we will use this structural property to manipulate the locomotor dynamics.

## 5. Manoeuvring by modulating sensitivity

In the previous sections, we showed how our CPG-based model can achieve a stable locomotion pattern, with the mechanical feedback term playing a key role in the emergence of travelling waves. The strength of this term is measured by the sensitivity strength $\sigma$, which was observed to influence the gait aspect as well as the direction of motion depending on its sign. Hence, time variations of the sensitivity function are likely to allow changes in locomotion, such as a direction reversal or turn.

The existence of interneurons controlling forward and backward locomotion in *C. elegans* has been experimentally well-established [60]. This command neurons are numerically studied and demonstrated to switch between forward and backward motions [61]. The turning motion of *C. elegans* is also well-known to be driven by interneural commands [62], although the detailed neuromechanical mechanisms are still under investigation.

In this section, we explore the behaviours obtained when following this idea and considering $\sigma$ as a function of time instead of a fixed parameter, taking inspiration from this *C. elegans* behaviour. The dynamical system governing the animal's motion now appears as a non-autonomous system, and its trajectory can be guided by tuning $\sigma$ as a control function.
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Transient crawling locomotion from an unstable saddle (UPO) to a stable limit cycle (SPO) with $N = 10$. (a) Snapshots of the shape in the $xy$ plane are superimposed on the trajectory of the end of the first link. The shape of the locomotor is plotted in different colours to present the time evolution and the black dot denotes the end of the first link in each snapshot. The object, initially located near the UPO ($\lambda = 1.27$) in the state space, exhibits an eventual convergence to the SPO ($\lambda = 0.97$) after transient behaviour. (b,c) Kymographs of the angle and phase dynamics, respectively. The model parameters are the same as in figure 6. See also electronic supplementary material, Video S6.</figcaption></figure>

## 5.1. Reverse motion

Sudden changes of the sensitivity value $\sigma$ with time were reported in a neuromechanical model of *C. elegans* crawling motion as a plausible cause for motion reversal and turns [61].

From the point of view of the dynamical system, equation (2.2), a simple sign change of $\sigma$ at some time $t_0$ modifies the dynamics landscape and in particular the position of the global SPO, which we know robustly emerges for both $\sigma$ and $-\sigma$. Thus, forward-backward motion reversal naturally occurs upon this sensitivity change.

When the coupling strength $C$ is sufficiently large, the locomotion immediately converges to the global SPO of the new dynamical system, and the crawler exhibits a quick reverse motion, which agrees well with *C. elegans* observations [37].

When the coupling weakens, the transient dynamics between the motion reversal exhibit more complicated behaviours, as shown in figure 10, where we present the transient behaviours of crawling locomotion with the same parameter set as in figure 9. We used the SPO of figure 9 as the initial configuration and changed the sign of $\sigma$ as $\sigma = 12 \mapsto \sigma = -12$ at time $t_0 = 0$. Snapshots of the shape gait are shown, with the horizontal axis indicating the time evolution. The arrow indicates the movement direction, and the crawler moves in the $-y$ direction until it turns to a backward motion moving in the $+y$ direction with a transient regime without phase synchrony (electronic supplementary material, video S7). As plotted in figure 10b,c the kymographs of the angle and phase dynamics also show complex transient dynamics until the angle reformulates the backward travelling wave.

## 5.2. Turning behaviour

So far, we have observed forward and backward motion, as well as the transition between both when the sensitivity strength $\sigma$ is switched between a set value and its opposite. However, this abrupt switch can be modulated to allow manoeuvring and reaction to changes in the environment. This would particularly enable foraging strategies such as run-and-tumble or any kind of taxis.

To illustrate this, we take the omega-turn behaviour of *C. elegans* as an example and propose a way to reproduce it in our CPG model. The turning strategy of *C. elegans* on agar gels is called an 'omega turn' because it temporarily breaks the symmetry of its gait by curving towards one side in a shape that roughly mimics the Greek letter $\Omega$ [37,56].

To allow these occasional symmetry breaks within the CPG dynamics, we add a term in equation (2.2) governing the evolution of $\phi$, which may reflect a neuronal command:

$$ \dot{\phi}_i = \omega_0 + \sum_j C \sin(\phi_j - \phi_i) + \sigma_1(t) \cos(\phi_i) S_i(t) + \sigma_2(t) K \cos(\phi_i - \psi_i), \quad (5.1) $$

where $K$ and $\psi_i$ are set to induce the local angle $\alpha_i$ to curve into some angle $\alpha_i^\Omega$ and the amplitude of $\sigma_2$ measures this forcing strength. Details on how to design $K$ and $\psi_i$ are provided in electronic supplementary material.

To clarify the role of the two activation functions in equation (5.1), the factor term of $\sigma_1$ sustains periodic locomotion as discussed previously, while the factor term in $\sigma_2$ drives each $\phi_i$ to converge to a constant value in a way that immobilizes the model animal in a given shape measured by target angles $\alpha_i^\Omega$. The relative amplitude of $\sigma_1$ and $\sigma_2$ determines which of these two behaviours dominate the instantaneous dynamics (electronic supplementary material, video S8).

With this set-up, when tuning $\sigma_1$ and $\sigma_2$ with respect to time, for example, following periodic functions, we may observe alternating forward or backward motion and also omega-turn-like behaviours effectively reorienting the animal's direction of motion between forward and backward sections of the trajectory (top diagram of figure 11a), where we used the same parameter set of the crawling locomotion as in previous figures (electronic supplementary material, video S9).

A family of functions for $\sigma_1$ and $\sigma_2$ that may yield an interesting range of behaviours is given by

$$ \left\{ \begin{array}{l} \sigma_1(t) = \sigma_1^0 \cos^{2n_1+1} \left( \frac{t}{T_1} + \chi \right), \\ \sigma_2(t) = \sigma_2^0 \sin^{2n_2+1} \left( \frac{t}{T_2} \right). \end{array} \right. \quad (5.2) $$

The integers $n_1$ and $n_2$ should be relatively small to obtain the distributed activation, while large values simulate impulse-like behaviour. This latter option is suitable for the omega-turn activation. The corresponding activation function is plotted in the
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Transient crawling behaviour when the sign of signal strength $\sigma$ switches. The SPO in figure 9 was used for the initial condition. We used the same parameter set except for $\sigma$, which was switched to $-\sigma$ at time $t_0 = 0$. (a) Time series of the snapshots. The arrow indicates the movement direction, and the crawler moves in the $-y$ direction until it turns to reverse to a backward motion in the $+y$ direction. (b,c) Kymographs of the angle and phase dynamics, respectively. See also electronic supplementary material, video S7.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Simulated turning and foraging behaviour. (a) Representation of the control functions $\sigma_1(t)$, $\sigma_2(t)$ in the $(\sigma_1, \sigma_2)$ plane. The top panel summarizes the different states preferred by the model animal at each extremity of the activation pattern: (clockwise from the top) omega turn, forward motion, mirrored omega turn and backward motion. The middle and bottom panels, respectively, display the activation functions used to generate the trajectories of panels (b) and (c). (b) Detailed aspect of the omega-turn simulation, with associated angle and phase kymographs. (c) Long-time behaviour with decorrelated frequencies for $\sigma_1$ and $\sigma_2$, showing complex, irregular trajectories. See also electronic supplementary material, Videos S8-S9.</figcaption></figure>

$\sigma_1 - \sigma_2$ plane in the middle panel of figure 11a, and the resulting short-term trajectory is shown in figure 11b. The snapshots of the shape gait are illustrated in different colours from blue to green to red as time progresses. The model animal initially moves towards the right, but then briefly turns and then starts moving again towards the left.

This mechanism can be exploited to simulate the animal exploring its environment. Further variability may be obtained by desynchronizing the periodicity of $\sigma_1$ and $\sigma_2$, taking, for example, $T_2/T_1 = \sqrt{2}$ in equation (5.2); the corresponding activation function is shown in the bottom plot of figure 11a, and the resulting long-term trajectory is displayed in figure 11c. Sharp turns at seemingly unpredictable angles due to the combination of the omega-turn mechanism and complex activation functions makes the trajectory strongly reminiscent of foraging behaviour [37].

## 6. Concluding remarks

In this study, we formulated a minimal model of neuromechanical undulatory locomotion of an elastic slender body in a dissipation-dominated environment using the RFT and coupled phase oscillators that form a CPG-like network. In contrast to existing CPG models, which assume an undulatory body motion, our network model preserves the motion reciprocity. Nonetheless, local phase modulation of each oscillator gives rise to robust undulatory locomotion in a self-organized manner, through mechanosensory feedback of the local mechanical load.

By numerical explorations, we obtained a stable, coherent undulatory motion in a broad parameter region and the emergent gait pattern reasonably agrees with various motions of *C. elegans* in different media. This process emerges from local symmetry breaking in the fore-aft resistive torque perceived by mechanoreceptors, which is then reflected at the network level by symmetry breaking with respect to the swimmer's transverse plane, allowing non-reciprocal movement. This elementary mechanism thus explains the link between local and global dynamics at the origin of the motor activation wave, a property that could be extended to other sensory feedback, such as proprioception.

This coherent undulatory locomotion was further explored using Purcell's three-link model. We found through theoretical and numerical analyses that a saddle-node bifurcation of periodic orbits occurs at a certain level of phase coupling. This indicates that the local CPG oscillators require a certain level of coupling to stabilize the network and that the sensory feedback cannot generate a stable gait on its own.
We then exploited this dynamical system structure to design manoeuvring strategies, by using the signal sensitivity as a control. Motivated by *C. elegans* behaviour on agar, we demonstrated through numerical explorations that our CPG model reproduces well the reverse motion and foraging behaviour known as the omega turn. Our deterministic model can thus include various predefined stereotyped behaviours, the so-called neuromechanical 'templates' [63] given by cycles and fixed points in the state space. Similarly, our model exhibits apparently random motion due to transient dynamics between these attractors. Recent data-driven modelling approaches may discover such attractors and the associated transients in the behavioural-state space [36,64,65].

One can also build upon this property to investigate, with a set of simple building blocks, complex long-term behaviours observed in some biological locomotors. Furthermore, this methodology of manipulating the nonlinear dynamical systems might offer a robust design method for redundant robot controllers [66].

In this study, we propose a general mechanism for gait adaptation through exteroceptive mechanosensory feedback, and remarkably, our numerical results are qualitatively in agreement with *C. elegans* behaviours, in particular on the decrease of wavelength of the emergent undulatory motion and amplification of internal actuation in high viscosity medium [52]. The current findings will be useful in modelling complex behavioural patterns and gait adaptations of animals.

The similar tendency of the changes in wavelength and internal actuation at the high viscosity regime is known in sperm motility [54,55]. In unicellular organisms with cilia and flagella, the undulatory motions are regulated by dynein motors, and mechanochemical oscillations of dynein actuation are reasonably represented by a coupled-oscillator system [67,68]. The exteroceptive feedback used in our manuscript could be used as a coarse-graining description of dynein regulation, although further studies should be warranted.

In our CPG model, the mechanosensation leads to a phase modulation that breaks the kinematic reversibility in the dissipative system. The non-reciprocity of the shape gait is recently found as the emergence of odd elasticity [43,69]. Such non-reciprocal material can generate effective work, and its mechanical efficiency merits further study.

In conclusion, our dynamical systems viewpoint is useful not only for understanding the extent of complexity and diversity found in animal locomotion but also for designing self-propelled, self-adaptive robots for a wide range of dissipative environments.

**Ethics.** This work did not require ethical approval from a human subject or animal welfare committee.

**Data accessibility.** The original code is freely available from the HAL online repository [70].

Supplementary material is available online [71].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** K.I.: conceptualization, investigation, software, visualization, writing—original draft, writing—review and editing; C.M.: conceptualization, investigation, software, visualization, writing—original draft, writing—review and editing; J.H.: conceptualization, investigation, software, visualization, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** K.I. acknowledges the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant No. 21H05309 and 24K21517) and the Japan Science and Technology Agency (JST), FOREST (Grant No. JPMJFR212N). C.M. acknowledges funding support for JSPS International Research Fellow (PE22023) (Grant No. 22KF0197). J.H. was supported by the French National Research Agency ANR (Grant No. ANR-2019-CE33-0004-01).

**Acknowledgements.** K.I. and C.M. were supported in part by the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located at Kyoto University. The authors acknowledge Prof. Toshiyuki Nakagaki for fruitful discussion.

## References

1. Nosrati R, Driouchi A, Yip CM, Sinton D. 2015 Two-dimensional slither swimming of sperm within a micrometre of a surface. *Nat. Commun.* **6**, 8703. (doi:10.1038/ncomms9703)

2. Gazzola M, Argentina M, Mahadevan L. 2014 Scaling macroscopic aquatic locomotion. *Nat. Phys.* **10**, 758–761. (doi:10.1038/nphys3078)

3. Jayne BC. 2020 What defines different modes of snake locomotion? *Integr. Comp. Biol.* **6**, 156–170. (doi:10.1093/icb/icaa017)

4. Di Santo V, Goerig E, Wainwright DK, Akanyeti O, Liao JC, Castro-Santos T, Lauder GV. 2021 Convergence of undulatory swimming kinematics across a diversity of fishes. *Proc. Natl Acad. Sci. USA* **118**, e2113206118. (doi:10.1073/pnas.2113206118)

5. Guasto JS, Estrada JB, Menolascina F, Burton LJ, Patel M, Franck C, Hosoi AE, Zimmer RK, Stocker R. 2020 Flagellar kinematics reveals the role of environment in shaping sperm motility. *J. R. Soc. Interface* **17**, 20200525. (doi:10.1098/rsif.2020.0525)

6. Gaffney EA, Ishimoto K, Walker BJ. 2021 Modelling motility: the mathematics of spermatozoa. *Front. Cell Dev. Biol.* **9**, 710825. (doi:10.3389/fcell.2021.710825)

7. Ijspeert AJ. 2008 Central pattern generators for locomotion control in animals and robots: a review. *Neural Networks* **21**, 642–653. (doi:10.1016/j.neunet.2008.03.014)

8. Wen Q, Gao S, Zhen M. 2018 *Caenorhabditis elegans* excitatory ventral cord motor neurons derive rhythm for body undulation. *Phil. Trans. R. Soc. B373*, 20170370. (doi:10.1098/rstb.2017.0370)

9. Thandiackal R et al. 2021 Emergence of robust self-organized undulatory swimming based on local hydrodynamic force sensing. *Sci. Robot.* **6**, f6354. (doi:10.1126/scirobotics.abf6354)

10. Tsybina YuA, Gordleeva SYu, Zharinov AI, Kastalski IA, Ermolaeva AV, Hramov AE, Kazantsev VB. 2022 Toward biomorphic robotics: a review on swimming central pattern generators. *Chaos Solitons Fractals* **165**, 112864. (doi:10.1016/j.chaos.2022.112864)

11. Wyart C, Carbo-Tano M. 2023 Design of mechanosensory feedback during undulatory locomotion to enhance speed and stability. *Curr. Opin. Neurobiol.* **83**, 102777. (doi:10.1016/j.conb.2023.102777)

12. Tytell ED, Hsu CY, Williams TL, Cohen AH, Fauci LJ. 2010 Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming. *Proc. Natl Acad. Sci. USA* **107**, 19832–19837. (doi:10.1073/pnas.1101564107)
13. Hamlet C, Fauci L, Morgan JR, Tytell ED. 2023 Proprioceptive feedback amplification restores effective locomotion in a neuromechanical model of lampreys with spinal injuries. *Proc. Natl Acad. Sci. USA* **120**, e2213302120. (doi:10.1073/pnas.2213302120)

14. Gordleeva SY, Kastalski JA, Tsybina YA, Yermolaeva AV, Hramov AE, Kazantsev VB. 2023 Control of movement of underwater swimmers: animals, simulated animates and swimming robots. *Phys. Life Rev.* **47**, 211–244. (doi:10.1016/j.plrev.2023.10.037)

15. Buchanan JT, Cohen AH. 1982 Activities of identified interneurons, motoneurons, and muscle fibers during fictive swimming in the lamprey and effects of reticulospinal and dorsal cell stimulation. *J. Neurophysiol.* **47**, 948–960. (doi:10.1152/jn.1982.47.5.948)

16. Di Prisco GV, Pearlstein E, Le Ray D, Robitaille R, Dubuc R. 2000 A cellular mechanism for the transformation of a sensory input into a motor command. *J. Neurosci.* **20**, 8169–8176. (doi:10.1523/jneurosci.20-21-08169.2000)

17. Ristroph L, Liao JC, Zhang J. 2015 Lateral line layout correlates with the differential hydrodynamic pressure on swimming fish. *Phys. Rev. Lett.* **114**, 018102. (doi:10.1103/PhysRevLett.114.018102)

18. Ishimoto K, Gadêlha H, Gaffney EA, Smith DJ, Kirkman-Brown J. 2017 Coarse-graining the fluid flow around a human sperm. *Phys. Rev. Lett.* **118**, 124501. (doi:10.1103/physrevlett.118.124501)

19. Hosoi AE, Goldman DI. 2015 Beneath our feet: strategies for locomotion in granular media. *Annu. Rev. Fluid Mech.* **47**, 431–453. (doi:10.1146/annurev-fluid-010313-141324)

20. Rieser JM, Gong C, Astley HC, Schiebel PE, Diaz K, Chong B, Pierce C *et al*. 2019 Geometric phase predicts locomotion performance in undulating living systems across scales. *arXiv* 1906.11374. (doi:10.48550/arXiv.1906.11374)

21. Florida Museum of Natural History. Biodiversity occurrence data. See https://www.invertebase.org/portal/taxa/index.php?taxauthid=1&taxon=86639&cidid=10.

22. Ding Y, Sharpe SS, Wiesenfeld K, Goldman DI. 2013 Emergence of the advancing neuromechanical phase in a resistive force dominated medium. *Proc. Natl Acad. Sci. USA* **110**, 10123–10128. (doi:10.1073/pnas.1302844110)

23. Treers LK, Cao C, Stuart HS. 2021 Granular resistive force theory implementation for three-dimensional trajectories. *IEEE Robot. Autom. Lett.* **6**, 1887–1894. (doi:10.1109/lra.2021.3057052)

24. Friedrich BM, Riedel-Kruse IH, Howard J, Jülicher F. 2010 High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory. *J. Exp. Biol.* **213**, 1226–1234. (doi:10.1242/jeb.039800)

25. Schulman RD, Backholm M, Ryu WS, Dalnoki-Veress K. 2014 Dynamic force patterns of an undulatory microswimmer. *Phys. Rev. E* **89**, 050701. (doi:10.1103/PhysRevE.89.050701)

26. Keaveny EE, Brown AEX. 2017 Predicting path from undulations for *C. elegans* using linear and nonlinear resistive force theory. *Phys. Biol.* **14**, 025001. (doi:10.1088/1478-3975/aa5ce6)

27. Dorgan KM, Law CJ, Rouse GW. 2013 Meandering worms: mechanics of undulatory burrowing in muds. *Proc. R. Soc. B* **280**, 20122948. (doi:10.1098/rspb.2012.2948)

28. Purcell EM. 1977 Life at low Reynolds number. *Am. J. Phys.* **45**, 3–11. (doi:10.1119/1.10903)

29. Ishimoto K, Yamada M. 2012 A coordinate-based proof of the scallop theorem. *SIAM J. Appl. Math.* **72**, 1686–1694. (doi:10.1137/110853297)

30. Fouad AD *et al*. 2018 Distributed rhythm generators underlie *Caenorhabditis elegans* forward locomotion. *elife* **7**, e29913. (doi:10.7554/elife.29913)

31. Ji H, Fouad AD, Teng S, Liu A, Alvarez-Illera P, Yao B, Li Z, Fang-Yen C. 2021 Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in *Caenorhabditis elegans*. *elife* **10**, e69905. (doi:10.7554/eLife.69905)

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

33. Cohen N, Denham JE. 2019 Whole animal modeling: piecing together nematode locomotion. *Curr. Opin. Syst. Biol.* **13**, 150–160. (doi:10.1016/j.coisb.2018.12.002)

34. Denham JE, Ranner T, Cohen N. 2018 Signatures of proprioceptive control in *Caenorhabditis elegans* locomotion. *Phil. Trans. R. Soc. B* **373**, 20180208. (doi:10.1098/rstb.2018.0208)

35. Johnson CL, Lewis TJ, Guy R. 2021 Neuromechanical mechanisms of gait adaptation in *C. elegans*: relative roles of neural and mechanical coupling. *SIAM J. Appl. Dyn. Syst.* **20**, 1022–1052. (doi:10.1137/20m1346122)

36. Ahamed T, Costa AC, Stephens GJ. 2021 Capturing the continuous complexity of behaviour in *Caenorhabditis elegans*. *Nat. Phys.* **17**, 275–283. (doi:10.1038/s41567-020-01036-8)

37. Stephens GJ, Johnson-Kerner B, Bialek W, Ryu WS. 2008 Dimensionality and dynamics in the behavior of *C. elegans*. *PLoS Comput. Biol.* **4**, e1000028. (doi:10.1371/journal.pcbi.1000028)

38. Ishimoto K, Moreau C, Herault J. Robust undulatory locomotion via neuromechanical adjustments in a dissipative medium. *arXiv* 2405.01802. (doi:10.48550/arXiv.2405.01802)

39. Moreau C, Giraldi L, Gadêlha H. 2018 The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella. *J. R. Soc. Interface* **15**, 20180235. (doi:10.1098/rsif.2018.0235)

40. Walker BJ, Ishimoto K, Gadêlha H, Gaffney EA. 2019 Filament mechanics in a half-space via regularised Stokeslet segments. *J. Fluid Mech.* **879**, 808–833. (doi:10.1017/jfm.2019.723)

41. Doi M. 2013 *Soft matter physics*. Oxford, UK: Oxford Academic. (doi:10.1093/acprof:oso/9780199652952.001.0001)

42. Ishimoto K, Moreau C, Yasuda K. 2022 Self-organized swimming with odd elasticity. *Phys. Rev. E* **105**, 064603. (doi:10.1103/PhysRevE.105.064603)

43. Ishimoto K, Moreau C, Yasuda K. 2023 Odd elastohydrodynamics: non-reciprocal living material in a viscous fluid. *PRX Life* **1**, 023002. (doi:10.1103/prxlife.1.023002)

44. Cohen AH, Ermentrout GB, Kiemel T, Kopell N, Sigvardt KA, Williams TL. 1992 Modelling of intersegmental coordination in the lamprey central pattern generator for locomotion. *Trends Neurosci.* **15**, 434–438. (doi:10.1016/0166-2236(92)90006-t)

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

46. Kuramoto Y. 1984 *Chemical oscillations, waves, and turbulence*. Berlin, Germany: Springer.

47. Kobayashi R, Tero A, Nakagaki T. 2006 Mathematical model for rhythmic protoplasmic movement in the true slime mold. *J. Math. Biol.* **53**, 273–286. (doi:10.1007/s00285-006-0007-0)

48. Kobayashi R, Nakagaki T, Ishiguro A. 2013 A design principle of the decentralized control and its applications. *AIP Conf. Proc.* **1558**, 2440–2443. (doi:10.1063/1.4826034)

49. Lighthill J. 1976 Flagellar hydrodynamics. *SIAM Rev.* **18**, 161–230. (doi:10.1137/1018040)

50. Gadêlha H, Gaffney EA, Smith DJ, Kirkman-Brown JC. 2010 Nonlinear instability in flagellar dynamics: a novel modulation mechanism in sperm migration? *J. R. Soc. Interface* **7**, 1689–1697. (doi:10.1098/rsif.2010.0136)

51. Sznitman J, Purohit PK, Krajacic P, Lamitina T, Arratia PE. 2010 Material properties of *Caenorhabditis elegans* swimming at low Reynolds number. *Biophys. J.* **98**, 617–626. (doi:10.1016/j.bpj.2009.11.010)

52. Fang-Yen C, Wyart M, Xie J, Kawai R, Kodger T, Chen S, Wen Q, Samuel ADT. 2010 Biomechanical analysis of gait adaptation in the nematode *Caenorhabditis elegans*. *Proc. Natl Acad. Sci. USA* **107**, 20323–20328. (doi:10.1073/pnas.1003016107)
53. Shen XN, Sznitman J, Krajacic P, Lamitina T, Arratia PE. 2012 Undulatory locomotion of *Caenorhabditis elegans* on wet surfaces. *Biophys. J.* **102**, 2772–2781. (doi:10.1016/j.bpj. 2012.05.012)

54. Smith DJ, Gaffney EA, Gadêlha H, Kapur N, Kirkman-Brown JC. 2009 Bend propagation in the flagella of migrating human sperm, and its modulation by viscosity. *Cell Motil. Cytoskelet.** **66**, 220–236. (doi:10.1002/cm.20345)

55. Ishimoto K, Gadêlha H, Gaffney EA, Smith DJ, Kirkman-Brown J. 2018 Human sperm swimming in a high viscosity mucus analogue. *J. Theor. Biol.* **446**, 1–10. (doi:10.1016/j.jtbi. 2018.02.013)

56. Pierce-Shimomura JT, Chen BL, Mun JJ, Ho R, Sarkis R, McIntire SL. 2008 Genetic analysis of crawling and swimming locomotory patterns in *C. elegans*. *Proc. Natl Acad. Sci. USA* **105**, 20982–20987. (doi:10.1073/pnas.0810359105)

57. Saiki Y. 2007 Numerical detection of unstable periodic orbits in continuous-time dynamical systems with chaotic behaviors. *Nonlinear Process. Geophys.* **14**, 615–620. (doi:10.5194/ npg-14-615-2007)

58. Ishimoto K, Gaffney EA. 2014 A study of spermatozoan swimming stability near a surface. *J. Theor. Biol.* **360**, 187–199. (doi:10.1016/j.jtbi.2014.06.034)

59. Hatton RL, Ding Y, Choset H, Goldman DI. 2013 Geometric visualization of self-propulsion in a complex medium. *Phys. Rev. Lett.* **110**, 078101. (doi:10.1103/PhysRevLett.110. 078101)

60. Chalfie M, Sulston J, White J, Southgate E, Thomson J, Brenner S. 1985 The neural circuit for touch sensitivity in *Caenorhabditis elegans*. *J. Neurosci.* **5**, 956–964. (doi:10.1523/ jneurosci.05-04-00956.1985)

61. Olivares EO, Izquierdo EJ, Beer RD. 2018 Potential role of a ventral nerve cord central pattern generator in forward and backward locomotion in *Caenorhabditis elegans*. *Netw. Neurosci.* **2**, 323–343. (doi:10.1162/netn_a_00036)

62. Kato S, Kaplan HS, Schrödel T, Skora S, Lindsay TH, Yemini E, Lockery S, Zimmer M. 2015 Global brain dynamics embed the motor command sequence of *Caenorhabditis elegans*. *Cell* **163**, 656–669. (doi:10.1016/j.cell.2015.09.034)

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

64. Gilpin W. 2020 Deep reconstruction of strange attractors from time series. *Adv. Neural Inf. Process. Syst.* **33**, 204–216.

65. Valani RN, Paganin DM. 2023 Attractor-driven matter. *Chaos* **33**, 023125. (doi:10.1063/5.0107401)

66. Seeja G, Selvakumar Arockia Doss A, Hency VB. 2022 A survey on snake robot locomotion. *IEEE Access* **10**, 112100–112116. (doi:10.1109/access.2022.3215162)

67. Cass JF, Bloomfield-Gadêlha H. 2023 The reaction-diffusion basis of animated patterns in eukaryotic flagella. *Nat. Commun.* **14**, 5638. (doi:10.1038/s41467-023-40338-2)

68. Maggi C, Saglimbeni F, Sosa VC, Di Leonardo R, Nath B, Puglisi A. 2023 Thermodynamic Limits of Sperm Swimming Precision. *PRX Life* **1**, 013003. (doi:10.1103/prxlife.1.013003)

69. Fruchart M, Scheibner C, Vitelli V. 2023 Odd viscosity and odd elasticity. *Annu. Rev. Condens. Matter Phys.* **14**, 471–510. (doi:10.1146/annurev-conmathphys-040821-125506)

70. Moreau C, Herault J, Ishimoto K. 2024 Code for the article: Robust undulatory locomotion via neuromechanical adjustments in a dissipative medium. *HAL*. See https://hal.science/ hal-04778432.

71. Ishimoto K, Moreau C, Herault J. 2025 Supplementary material from: Robust undulatory locomotion via neuromechanical adjustments in a dissipative medium. Figshare. (doi:10. 6084/m9 figshare.c.7639319)