Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in Caenorhabditis elegans

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

Neural circuits coordinate with muscles and sensory feedback to generate motor behaviors appropriate to an animal’s environment. In C. elegans, the mechanisms by which the motor circuit generates undulations and modulates them based on the environment are largely unclear. We quantitatively analyzed C. elegans locomotion during free movement and during transient optogenetic muscle inhibition. Undulatory movements were highly asymmetrical with respect to the duration of bending and unbending during each cycle. Phase response curves induced by brief optogenetic inhibition of head muscles showed gradual increases and rapid decreases as a function of phase at which the perturbation was applied. A relaxation oscillator model based on proprioceptive thresholds that switch the active muscle moment was developed and is shown to quantitatively agree with data from free movement, phase responses, and previous results for gait adaptation to mechanical loadings. Our results suggest a neuromuscular mechanism underlying C. elegans motor pattern generation within a compact circuit.

Introduction

Animal display locomotor behaviors such as crawling, walking, swimming, or flying via rhythmic patterns of muscle contractions and relaxations. In many animals, motor rhythms originate from networks of central pattern generators (CPGs), neuronal circuits capable of generating rhythmic outputs without rhythmic input (Cohen and Wallen, 1980; Grillner, 2003; Kiehn, 2011; Kristan and Calabrese, 1976; Marder and Calabrese, 1996; Pearce and Friesen, 1984; Yu et al., 1999). CPGs typically generate rhythms through reciprocal inhibitory synaptic interactions between two populations. In vertebrates, motor rhythms arise from half-center oscillator modules in the spinal cord (Marder and Calabrese, 1996).

Although isolated CPGs can produce outputs in the absence of sensory input, in the intact animal sensory feedback plays a critical role in coordinating motor rhythms across the body and modulating their characteristics (Friesen, 2009; Grillner and Wallén, 2002; Mullins et al., 2011; Pearson, 2004; Wen et al., 2012). Sensory feedback allows animals to adapt locomotor patterns to their surroundings (Andersson et al., 1981; Brodfuehrer and Friesen, 1986) and adapt to unexpected perturbations (Ekeberg and Grillner, 1999). In leeches (Cang et al., 2001; Cang and Friesen, 2000) and Drosophila (Akitake et al., 2015; Mendes et al., 2013), specialized proprioceptive neurons and sensory receptors in body muscles detect sensory inputs to regulate and coordinate the centrally generated motor patterns. In limbed vertebrates, proprioceptors located in muscles, joints, and/or skin detect body movements and interact with premotor interneurons to coordinate limb movements (Pearson, 2004). Sensory inputs induced by electric stimulation of receptor cells (Yu and Friesen, 2004) or by mechanical perturbation of body segments (Grillner et al., 1981) can entrain an animal’s motor behavior to imposed patterns, demonstrating the flexibility of motor systems in responding to feedback.

Animal movements are driven not only by active muscle contractions but also by passive mechanical forces including elastic recoil of muscles and other body structures, internal damping forces, and forces from the interaction with the external environment. Efficient locomotion in vertebrates depends on storage of elastic energy in tendons and muscles (Roberts and Azizi, 2011). In insects, elasticity in the leg joint plays an important role in generating forces for walking and jumping (Ache and Matheson, 2013). A comprehensive understanding of animal locomotion should therefore encompass not only neural activity, muscle activity, and sensory feedback, but also biomechanical forces within the animal’s body and between the animal and its environment (Figure 1A; Borgmann et al., 2009; Grillner and Wallén, 2002; Kiehn, 1998).

Figure 1

Download asset Open asset

Rhythm generation in *C. elegans*.

(A) Motor neurons generate neuronal signals to control the activation of body wall muscles (BWM), which generates movement subject to internal and external environmental constraints. Sensory input provides feedback about body position and the environment. (**B,C**) Two possible models for locomotory rhythm generation in *C. elegans*. (B) In a reflex loop model, sensory neurons (SN) detect body postures and excite motor neurons (MN) to activate body wall muscles. (C) In a central pattern generator (CPG) model, network of motor neurons generates basic rhythmic patterns that are transmitted to body wall muscles while sensory feedback modulates the CPG rhythm. Diagrams (**B-C**) are adapted from Figure 1 in Marder and Bucher, 2001.

Here, we study mechanisms of locomotor rhythm generation and its modulation by sensory feedback in the nematode Caenorhabditis elegans. With its easily quantifiable behavior (Croll, 1971), well-mapped nervous system (Cook et al., 2019; White et al., 1986), genetic manipulability (Bargmann, 1998; Brenner, 1974; Hobert, 2003), and optical transparency, this worm is a unique model for obtaining an integrative understanding of locomotion.

C. elegans forward locomotion consists of anterior-to-posterior dorsoventral undulations (Croll, 1971). These movements are mediated by a neuromuscular circuit consisting of interneurons, excitatory cholinergic motor neurons, inhibitory GABAergic motor neurons, and body wall muscles. Laser ablation studies have shown that the cholinergic B-type motor neurons are required for forward locomotion (Chalfie et al., 1985). The GABAergic D-type motor neurons provide dorsoventral cross-inhibition to the body wall muscles and are essential for maintaining normal wave shape and frequency during forward locomotion (Deng et al., 2021; McIntire et al., 1993). A set of premotor interneurons (AVB, PVC, AVA, AVD, and AVE) regulate forward and reverse movements (Chalfie et al., 1988; Driscoll and Kaplan, 1997; Von Stetina et al., 2006). Ablation of all premotor interneurons or the D-type motor neurons does not deprive C. elegans of the ability to move forward (Chalfie et al., 1985; Gao et al., 2018; Kawano et al., 2011), suggesting that a network consisting of excitatory motor neurons and muscles may be sufficient to generate rhythmicity. Optogenetic and lesion experiments have shown that multiple oscillators exist in the ventral nerve cord (Fouad et al., 2018). However, the mechanisms that give rise to these oscillators are still poorly understood.

Proprioceptive feedback is crucial for C. elegans motor behavior. Studies have identified several neuron classes that have proprioceptive roles. The B-type motor neurons mediate proprioceptive coupling between anterior to posterior bending during forward locomotion (Wen et al., 2012). The SMDD motor neurons, localized at the head, have been identified as proprioceptive regulators of head steering during locomotion (Yeon et al., 2018). Both the B-type motor neurons and the SMDD head motor neurons have long asynaptic processes hypothesized to have proprioceptive function (White et al., 1986) and have been suggested as candidate locomotor CPG elements (Kaplan et al., 2020). In addition, two types of neurons, the DVA and PVD interneurons, have also been described as having proprioceptive roles in the regulation of worm’s body bend movement. The cell DVA has been shown to exhibit proprioceptive properties with a dependence on a mechanosensitive channel, TRP-4, which acts as a stretch receptor to regulate the body bend amplitude during locomotion (Li et al., 2006). In another study, body bending was shown to induce local dendritic calcium transients in PVD and dendritic release of a neuropeptide encoded by nlp-12, which appears to regulate the amplitude of body movements (Tao et al., 2019).

To experimentally probe mechanisms of rhythmic motor generation, including the role of proprioceptive feedback, we measured the phase response curve (PRC) upon transient optogenetic inhibition of the head muscles. We found that the worms displayed a biphasic, sawtooth-shaped PRC with sharp transitions from phase delay to advance.

We used these findings to develop a computational model of rhythm generation in the C. elegans motor circuit in which a relaxation-oscillation process, with switching based on proprioceptive feedback, underlies the worm’s rhythmic dorsal-ventral alternation. Computational models for C. elegans motor behavior have long been an important complement to experimental approaches, since an integrative understanding of locomotion requires consideration of neural, muscular, and mechanical degrees of freedom, and are often tractable only by modeling (Boyle et al., 2012; Bryden and Cohen, 2008; Denham et al., 2018; Izquierdo and Beer, 2018; Johnson et al., 2021; Karbowski et al., 2008; Kunert et al., 2017; Olivares et al., 2021). We sought to develop a phenomenological model to describe an overall mechanism of rhythm generation but not the detailed dynamics of specific circuit elements. We aimed to incorporate biomechanical constraints of the worm’s body and its environment (Fang-Yen et al., 2010; Gray and Lissmann, 1964; Wallace, 1968), as well as account for how sensory feedback is incorporated. To improve predictive power, we aimed to minimize the number of free parameters used in the model. Finally, we sought to optimize and test this model with new experiments as well as with published findings.

Our model reproduces the observed PRC and describes the locomotory dynamics around optogenetic inhibitions in a manner that closely fits our experimental observations. Our model also agrees with results on gait adaptation to external load and the asymmetry in time-dependent curvature patterns of undulating worms. Our experimental findings and computational model together yield insights into how C. elegans generates rhythmic locomotion and modulates them depending on the environment.

Results

C. elegans forward locomotion exhibits a stable and nonsinusoidal limit cycle

To gain insight into wave generation, we first sought to examine the quantitative behavioral characteristics of worms during forward locomotion. First, we measured the undulatory dynamics of body bending by computing the time-varying curvature along the centerline of the body (Fang-Yen et al., 2010; Leifer et al., 2011; Pierce-Shimomura et al., 2008; Wen et al., 2012) from analysis of dark field image sequences of worms exhibiting forward locomotion. In order to quantitatively treat the drag between the body and its environment, we examined locomotion of worms in dextran solutions of known viscosity (see Appendix; Fang-Yen et al., 2010). The normalized body coordinate is defined by the distance along the body centerline divided by the body length (Figure 2A). The curvature $κ$ at each point along the centerline of the body is the reciprocal of local radius of curvature (Figure 2A), with a positive (negative) curvature representing ventral (dorsal) bending. We further define the dimensionless or scaled curvature $K = κ \cdot L$ , where $L$ is the length of the worm. Using this metric, we quantified the worm’s forward movement by calculating scaled curvature as a function of body coordinate and time (Figure 2B).

Figure 2 with 1 supplement see all

Download asset Open asset

Undulatory dynamics of freely moving worms.

(A) Worm undulatory dynamics are quantified by the time-varying curvature along the body. The normalized body coordinate is defined by the fractional distance along the centerline (head = 0, tail = 1). The curvature $κ$ is the reciprocal of the local radius of curvature with positive and negative values representing dorsal and ventral curvature, respectively. (B) Curvature as a function of time and body coordinate during forward movement in a viscous liquid. Body bending curvature $K$ is represented using the nondimensional product of $κ$ and body length $L$ . (C) Curvature (black) in the anterior region (average over body coordinate 0.1-0.3) and the time derivative of curvature (dashed purple). Red circles mark four representative phases (0, $π / 2$ , $π$ , and $3 π / 2$ ). The curve is an average of 5041 locomotory cycles from 116 worms. (D) Phase portrait representation of the oscillatory dynamics of the anterior region, showing the curvature and the time derivative of the curvature parameterized by time. Images of worm correspond to the phases marked in (C). Arrow indicates clockwise movement over time. Dash-boxed region of the worm body indicates the 0.1–0.3 body coordinates. h: head; t: tail; v: ventral side; d: dorsal side. Gray curves are individual locomotory cycles from freely moving worms (10 randomly selected cycles are shown). (*Inset*) waveform of the scaled active muscle moment, estimated by equation $M_{a} = K + τ_{u} \dot{K}$ . Both curves were computed from the data used in (C).

We used this behavioral data to generate phase portraits, geometric representations of a dynamical system’s trajectories over time (Izhikevich, 2007), in which the time derivative of the curvature is plotted against the curvature. If the curvature were sinusoidal over time, as it is often modeled in slender swimmers (Fang-Yen et al., 2010; Gray, 1933; Guo and Mahadevan, 2008; Niebur and Erdös, 1991), the time derivative of curvature would also be sinusoidal, with a phase shift of $π / 4$ radians relative to the curvature, and the resulting phase portrait would be symmetric about both the $K$ and $d K / d t$ axes. Instead, we found that the phase portrait of the bending movement in the worm’s head region (0.1–0.3 body coordinate) during forward locomotion is in fact non-ellipsoidal and strongly asymmetric with respect to reflection across the $K$ or $d K / d t$ axes (Figure 2D). Plots of both the phase portrait (Figure 2D) and the time dependence (Figure 2C) show that $K$ and $d K / d t$ are strongly non-sinusoidal.

In addition to the head, other parts of the worm’s body also display nonsinusoidal bending movements (Figure 2—figure supplement 1). In this paper, we focus on curvature dynamics of the worm’s head region (0.1–0.3 body coordinate) where the bending amplitude is largest and the nonsinusoidal features are most prominent (Figure 2—figure supplement 1).

We asked whether the phase portrait represents a stable cycle, that is whether the system tends to return to the cycle after fluctuations or perturbations away from it. To this end, we analyzed the recovery after brief optogenetic muscle inhibition. We used a closed-loop system for optically targeting specific parts of the worm (Fouad et al., 2018; Leifer et al., 2011) to apply brief pulses of laser illumination (0.1 s duration, 532 nm wavelength) to the heads of worms expressing the inhibitory opsin NpHR in body wall muscles (via the transgene Pmyo-3::NpHR). Simultaneous muscle inhibition on both sides causes C. elegans to straighten due to internal elastic forces (Fang-Yen et al., 2010). Brief inhibition of the head muscles during forward locomotion was followed by a maximum degree of paralysis approximately 0.3 s after the end of the pulse, then a resumption of undulation (Figure 3A,B; Video 1).

Figure 3 with 8 supplements see all

Download asset Open asset

Analysis of phase-dependent inhibitions for head oscillation using transient optogenetic muscle inhibition.

(A) Images of a transgenic worm (*Pmyo-3::NpHR*) perturbed by a transient optogenetic muscle inhibition in the head during forward locomotion. Green shaded region indicates the 0.1 s laser illumination interval. h: head; t: tail; v: ventral side; d: dorsal side. (B) Effect of muscle inhibition on mean absolute curvature of the head. Black curve represents control ATR+ (no light) group (3523 measurements using 337 worms). Brown curve represents control ATR- group (2072 measurements using 116 worms). Red curve represents ATR+ group (1910 measurements using 337 worms). Green bar indicates 0.1 s light illumination interval starting at $t = 0$ . (**C-E**) Perturbed dynamics around light pulses occurring in the phase range $[0, π / 6]$ . (C) Kymogram of time-varying curvature $K$ around a 0.1 s inhibition (green dashed box). (D) Mean curvature dynamics around the inhibitions (green bar, aligned at $t = 0$ ) from ATR+ group (red curve, 11 trials using 4 worms) and control ATR+ (no light) group (black curve, eight trials using three worms). Gray curves are individual trials from ATR+ group (10 randomly selected trials are shown). (E) Mean phase portrait graphs around the inhibitions (green line) from ATR+ group (same trials as in D) and control group (ATR+, no light, 3998 trials using 337 worms). Gray curves are individual trials from ATR+ group. (**F-H**) Similar to (**C-E**), for phase range $[π / 3, π / 2]$ . (**I-K**) Similar to (**C-E**), for phase range $[2 π / 3, 5 π / 6]$ . (L) PRC from optogenetic inhibition experiments (ATR+ group, 991 trials using 337 worms, each point indicating a single illumination of one worm). The curve was obtained via a moving average along the x-axis with $0.16 π$ in bin width and the filled area represents 95% confidence interval within the bin. (M) A 2-dimensional histogram representation of the PRC using the same data. The histogram uses 25 bins for both dimensions, and the color indicates the number of data points within each rectangular bin.

Video 1

Download asset

posterframe for video — Transient illumination of the anterior region of a freely moving *Pmyo-3::NpHR* worm.

Green-shaded region indicates timing and location of illumination.

To quantify the recovery dynamics, we defined a normalized deviation $d$ describing the state of the system relative to the phase portrait of normal oscillation (see Appendix), such that $d = - 1$ at the origin, $d = 0$ at the limit cycle, and $d > 0$ outside the limit cycle. We found that the deviation following optogenetic perturbation (Figure 3—figure supplement 1) decays toward zero regardless of the initial deviation from the normal cycle, indicating that the worm tends to return to its normal oscillation after a perturbation. These results show that C. elegans head oscillation during forward locomotion is stable under optogenetic perturbation. The dynamics of these perturbed worms also allow us to reconstruct the phase isochrons and vector flow fields (Figure 3—figure supplement 2) of the worm’s head oscillation, two other important aspects of an oscillator (see Appendix).

Taken together, these results show that during forward locomotion, head oscillation of a worm constitutes a stable oscillator containing a nonsinusoidal limit cycle.

Transient optogenetic inhibition of head muscles yields a slowly rising, rapidly falling phase response curve

The phase response curve (PRC) describes the change in phase of an oscillation induced by a perturbation as a function of the phase at which the perturbation is applied, and is often used to characterize biological and nonbiological oscillators (Izhikevich, 2007; Pietras and Daffertshofer, 2019; Schultheiss et al., 2011). We performed a phase response analysis of the worm’s locomotion upon transient optogenetic inhibitions.

Using data from 991 illuminations (each 0.1 s in duration) in 337 worms, we analyzed the animals’ recovery from transient paralysis as a function of the phase at which the illumination occurred. We define the phase such that it equals to zero at the point of maximum ventral bending (Figure 3D). When inhibition occurred with phase in the interval $[0, π / 6]$ , the head typically straightened briefly and then continued the previous bend, resulting in a phase delay for the oscillation (Figure 3C–E). When inhibition occurred with phase in the interval $[π / 3, π / 2]$ , the head usually appeared to discontinue the previous bend movement, which resulted in a small phase advance (Figure 3F–H). When inhibition occurred with phase in the interval $[2 π / 3, 5 π / 6]$ , the head response was similar to that within the interval $[0, π / 6]$ , and also resulted in a phase delay (Figure 3I–K).

Combining the data from all phases of inhibition yielded a sawtooth-shaped PRC with two sharp transitions from phase delay to advance as well as two relatively slow ascending transitions from phase advance to delay (Figure 3L,M). In control worms, which do not express NpHR in the body wall muscles (see Materials and methods), the resulting PRC shows no significant phase shift over any phases of illumination (Figure 3—figure supplement 3). In worms perturbed with shorter pulses (0.055 s duration), we observed a similar sawtooth-shaped PRC (Figure 3—figure supplement 4).

In addition to phase response analyses with perturbations to the worm’s anterior region, we conducted similar analyses for the dynamics across the body by optogenetically inhibiting body wall muscles of other regions (Figure 3—figure supplement 5). We found that the sawtooth feature of PRC tends to decrease monotonically as the perturbation occurs further away from the head (Figure 3—figure supplement 5A,E,I).

Next, we asked whether the sharp downward transitions in the PRC represent a continuous decrease or instead result from averaging data from a bimodal distribution. When we plotted the distribution of the same data in a 2-D representation we found that the phase shifts display a piecewise, linear increasing dependence on the phase of inhibition with two abrupt jumps occurring at $ϕ \approx π / 3$ and $4 π / 3$ , respectively (Figure 3M). This result shows that the sharp decreasing transitions in PRC reflect bimodality in the data rather than continuous transitions.

In addition to examining PRCs induced by muscle inhibition, we also calculated PRCs with respect to inhibitions of cholinergic motor neurons. We performed similar experiments on transgenic worms in which the inhibitory opsin NpHR is expressed in either all cholinergic neurons (Punc-17::NpHR::ECFP) or B-type motor neurons (Pacr-5::Arch-mCherry). In both strains, we again observed sawtooth-shaped PRCs (Figure 3—figure supplements 6 and 7), with variations only in the magnitudes of phase shifts. These experiments show that the sawtooth-shaped feature of PRC is maintained for motor neuron inhibition, suggesting that the transient muscle and neuron inhibition interrupt the motor circuit dynamics in a similar manner.

The GABAergic D-type motor neurons provide a dorsoventral reciprocal inhibition of opposing muscles during locomotion. We asked whether the D-type motor neurons are required for the observed sawtooth shape of the PRC. We examined transgenic worms that express NpHR in the body wall muscles but have mutations unc-49(e407), a loss-of-function mutant of GABA_A receptor that is required by the D-type motor neurons (Bamber et al., 1999). After performing optogenetic inhibition experiments, we found that the PRC also displays a sawtooth feature (Figure 3—figure supplement 8). This result shows that D-type motor neurons are not necessary for the motor rhythm generator to show the sawtooth-shaped PRC.

Sawtooth-shaped PRCs are observed in a number of systems with oscillatory dynamics, including the van der Pol oscillator (Cestnik and Rosenblum, 2018), and may reflect a phase resetting property of an oscillator with respect to a perturbation (Izhikevich, 2007; Schultheiss et al., 2011). Further interpretation of the PRC results is given below.

Worm muscles display a rapid switch-like alternation during locomotion

As a first step in interpreting and modeling our findings, we estimated the patterns of muscle activity in freely moving worms, in part by drawing on previous biomechanical analyses of nematode movement (Fang-Yen et al., 2010; Gray and Lissmann, 1964; Wallace, 1968).

In mechanics, a moment is a measure of the ability of forces to produce bending about an axis. Body wall muscles create local dorsal or ventral bending by generating active moments across the body. In addition to the active moments from muscles, there are also passive moments generated by the worm’s internal viscoelasticity and by the forces due to the interaction of the worm with its external environment.

We estimated the output patterns of the active muscle moment that drives the head oscillations of freely moving worms immersed in viscous solutions. Following previous analyses of C. elegans locomotor biomechanics under similar external conditions (Fang-Yen et al., 2010), the scaled active muscle moment can be described as a linear combination of the curvature and the time derivative of the curvature (Equation 1; also see Methods and Appendix). We observed that in the phase portrait graph (Figure 2D), there are two nearly linear portions of the curve. We hypothesized that these linear portions correspond to two bouts during which the active muscle moment is nearly constant.

Using fits to the phase plot trajectory (see Materials and methods and Appendix) we estimated the waveform of the active muscle moment as a function of time (Figure 2D Inset). We found that the net active muscle moment alternates between two plateau regions during forward locomotion. From the slope of the steep portions on this curve, we estimated the time constant for transitions between active moments to be $τ_{m} \approx 100 m s$ . This time constant is much smaller than the duration of each muscle moment plateau period ( $\approx 0.5 s$ ), suggesting that the system undergoes rapid switches of muscle contractions between two saturation states.

A relaxation oscillator model explains nonsinusoidal dynamics

We reasoned that the rapid transitions of the active muscle moment might reflect a switching mechanism in the locomotory rhythm generation system. We hypothesized that the motor system generates locomotory rhythms by switching the active moment of the muscles based on proprioceptive thresholds.

To expand further upon these ideas, we developed a quantitative model of locomotory rhythm generation. We consider the worm as a viscoelastic rod where the scaled curvature K(t) varies according to:

\begin{matrix} K (t) + τ_{u} \frac{d K (t)}{d t} = M_{a} (t), \end{matrix}

where $τ_{u}$ describes the time scale of bending relaxation and $M_{a} (t)$ is the time-varying active muscle moment scaled by the bending modulus and the body length (see detailed derivations in Appendix). We note that in a stationary state ( $d K / d t = 0$ ), the curvature would be equal to the scaled active muscle moment. That is, the scaled active moment represents the static curvature that would result from a constant muscle moment.

We define a proprioceptive feedback variable $P$ as a linear combination of the current curvature value and the rate of change of curvature. In our model, once this variable reaches either of two thresholds $P_{t h}$ and $- P_{t h}$ (Figure 4D), the active muscle moment undergoes a change of sign (Figure 4E), causing the head to bend toward the opposite direction (Figure 4B).

Figure 4 with 1 supplement see all

Download asset Open asset

Free-running dynamics of a bidirectional relaxation oscillator model.

(A) Schematic diagram of the relaxation oscillator model. In this model, sensory neurons (SN) detect the total curvature of the body segment as well as the time derivative of the curvature. The linear combination of the two values, $P = K + b \dot{K}$ , is modeled as the proprioceptive signal which is transmitted to motor neurons (MN). The motor neurons alternatingly activate dorsal or ventral body wall muscles (BWM) based on a thresholding rule: (1) if $P < - P_{t h}$ , the ventral body wall muscles get activated and contract while the dorsal side of muscles relax; (2) if $P > P_{t h}$ , vice versa. Hence, locomotion rhythms are generated from this threshold-switch process. (B) Time-varying curvature $K$ of the model oscillator. The time axis is normalized with respect to oscillatory period (same for D, E, and F). (C) Phase portrait graph of the model oscillator. Proprioceptive threshold lines (gray dashed lines) intersect with the phase portrait graph at two switch points (red circles) at which the active moment of body wall muscles is switched. (D) Time-varying proprioceptive feedback $P$ received by the motor neurons. Horizontal lines denote the proprioceptive thresholds (gray dashed lines) that switch the active muscle moment at switch points (red circles, intersections between the proprioceptive feedback curve and the threshold lines). (E) Time-varying active muscle moment. Blue-dashed square wave denotes target moment ( $M_{t}$ ) that instantly switches directions at switch points. Black curve denotes the active muscle moment ( $M_{a}$ ) which follows the target moment in a delayed manner. (F) Time varying curvature in the worm’s head region from experiments (red, 5047 cycles using 116 worms) and model (black). Model curvature matches experimental curvature with an MSE ≈ 0.18. (*Inset*) Bar graph of $U$ (time period of bending toward the ventral or dorsal directions) and $D$ (time period of straightening toward a straight posture). Vertical bars are averages of fractions with respect to undulatory period $T_{0}$ of $U$ and $T$ (*** indicates p<0.0005 using Student’s t test).

Our model has 5 parameters: (1) $τ_{u}$ , the bending relaxation time scale, (2) $τ_{m}$ , the muscle switching time scale, (3) $M_{0}$ , the amplitude of the scaled active muscle moment, (4-5) $b$ and $P_{t h}$ , which determine the switch threshold. The first three parameters were directly estimated from our experimental results from freely moving worms (see Appendix). Parameters $b$ and $P_{t h}$ were obtained using a two-round fitting procedure by fitting the model first to the freely moving dynamics (first round) and then to the experimental phase response curve (second round) (see Appendix).

With this set of parameters, we calculated the model dynamics as represented by the phase portrait (Figure 4C) as well as curvature waveform in one cycle period (Figure 4F). We found that in both cases the model result agreed with our experimental observations. Our model captures the asymmetric phase portrait trajectory shape found from our experiments (Figure 2D). It also describes the asymmetry of head bending during locomotion: bending toward the ventral or dorsal directions occurs slower than straightening toward a straight posture during the locomotory cycle (Figure 4F Inset).

Considering the hypothesized mechanism under the biomechanical background (Equation 1), our model provides a simple explanation for the observed bending asymmetry during locomotion. According to the model, the active muscle moment is nearly constant during each period between transitions of the muscle moment. Biomechanical analysis under this condition predicts an approximately exponential decay in curvature, which gives rise to an asymmetric feature during each half period (Figure 4F).

Relaxation oscillator model reproduces responses to transient optogenetic inhibition

We performed simulations of optogenetic inhibitions in our model. To model the transient muscle paralysis, the muscle moment is modulated by a bell-shaped function of time (Figure 4—figure supplement 1; also see Appendix) such that, upon inhibition, it decays toward zero and then recovers to its normal value, consistent with our behavioral observations (Figure 3B).

From simulations with different sets of model parameters, we found that the model PRCs consistently exhibited the sawtooth shape found in experiments, although differing in height and timing of the downward transitions. In addition to the model parameters $τ_{u}$ , $M_{0}$ , and $τ_{m}$ that had been explicitly estimated from free-moving experiments, we performed a two-round fitting procedure (see Appendix) to determine the other parameters (including $b$ , $P_{t h}$ , and parameters for describing the optogenetically induced muscle inhibitions (see Figure 4—figure supplement 1)) to best fit the freely moving dynamics and the experimental PRC, respectively, with a minimum mean squared error (MSE) (Figures 4F and 5A; also see Appendix). For the parameters $b$ and $P_{t h}$ , the optimization estimated their values to be $b = 0.046 s$ and $P_{t h} = 2.33$ , as shown on the phase portraits (gray dashed lines in Figures 4C, 5B and D).

Figure 5 with 1 supplement see all

Download asset Open asset

Simulations of optogenetic inhibitions in the relaxation oscillator model.

(A) Phase response curves measured from experiments (blue, same as in Figure 3L) and model (orange). Model PRC matches experimental PRC with an MSE ≈ 0.12. (**B,C**) Simulated dynamics of locomotion showing inhibition-induced phase delays in the model oscillator. (B) Simulated phase portrait graphs around inhibition occurring at $π / 6$ phase of cycle for perturbed (red) and unperturbed (black) dynamics. Green bar indicates the phase during which the inhibition occurs. (C) Same dynamics as in (B), represented by time-varying curvatures. The time axis is normalized with respect to oscillatory period (same for E). (**D,E**) Simulated dynamics of locomotion showing inhibition-induced phase advances in the model oscillator. (D) Simulated phase portrait graphs around inhibition occurring at $π / 2$ phase of cycle for perturbed (red) and unperturbed (black) dynamics. (E) Same dynamics as in (D), represented by time-varying curvatures.

The threshold-switch mechanism model provides an explanation for the observed sawtooth-shaped PRC. By comparing model phase portrait graphs around inhibitions occurring at different phases (Figure 5B–E), we found that the phase shift depends on the relative position of the inhibition with respect to the switch points on the phase plane. (1) If the effect of the inhibition occurs before the system reaches its switch point (Figure 5B), the system will recover by continuing the previous bend and the next switch in the muscle moment will be postponed, thereby leading to a phase delay (Figure 5C). (2) As the inhibition progressively approaches the switch point, one would expect that the next switch in the muscle moment will also be progressively postponed; this explains the increasing portions of the PRC. (3) If the inhibition coincides with the switch point (Figure 5D), the muscle moment will be switched at this point and the system will recover by aborting the previous bend tendency, resulting in a small phase advance (Figure 5E). This switching behavior explains the two sharp downward transitions in the PRC.

Relaxation oscillator model predicts phase response curves for single-side muscle inhibition

As a further test of the model, we asked what PRCs would be produced with only the ventral or dorsal head muscles being transiently inhibited. In the model, the muscle activity is represented using the scaled active moment of muscles. We conducted model simulations (see Appendix) to predict the PRCs for transient inhibitions of muscles on the dorsal side (Figure 6A, Upper) and ventral side (Figure 6B, Upper), respectively.

Figure 6 with 2 supplements see all

Download asset Open asset

The model predicts phase response curves with respect to single-side muscle inhibitions.

(A) (*Upper*) a schematic indicating a transient inhibition of body wall muscles of the head on the dorsal side. (*Lower*) the corresponding PRC measured from experiments (blue, 576 trials using 242 worms) and model (orange). (B) (*Upper*) a schematic indicating a transient inhibition of body wall muscles of the head on the ventral side. (*Lower*) the corresponding PRC measured from experiments (blue, 373 trials using 176 worms) and model (orange). For the two experiments, each point indicates a single illumination (0.1 s duration, 532 nm wavelength) of one worm. Experimental curves were obtained using a moving average along the x-axis with $0.16 π$ in bin width. Filled area of each experimental curve represents 95% confidence interval with respect to each bin of data points.

To experimentally perform phase response analysis of single-side muscle inhibitions, we visually distinguished each worm’s dorsoventral orientation (via vulval location) and targeted light to either the ventral or dorsal side of the animal. Transiently illuminating (0.1 s duration) dorsal or ventral muscles in the head region of the transgenic worms (Pmyo-3::NpHR) induced a brief paralyzing effect when the segment was bending toward the illuminated side but did not induce a significant paralyzing effect when the segment was bending away from the illuminated side (Figure 6—figure supplement 1).

Combining the experimental data from all phases of dorsal-side or ventral-side inhibition yielded the corresponding PRCs (Figure 6A,B, respectively), from which we found that both PRCs show a peak in the phase range during which the bending side is illuminated but shows no significant phase shift in the other phase range. The experimental observations are qualitatively consistent with model predictions.

We found that the PRC of dorsal-side illumination shows a smaller paralytic response than that of ventral-side illumination. This discrepancy may be due to different degrees of paralysis achieved during ventral vs. dorsal illumination (Figure 6—figure supplement 1), possibly due to differences in levels of opsin expression and/or membrane localization. We therefore modulated the parameter for describing degree of paralysis when simulating the PRC of the dorsal-side illumination to qualitatively account for this discrepancy (see Appendix).

Our model is consistent with the dependence of wave amplitude and frequency on external load

C. elegans can swim in water and crawl on moist surfaces, exhibiting different undulatory gaits characterized by different frequency, amplitude, and wavelength (Figure 7A). Previous studies Berri et al., 2009; Fang-Yen et al., 2010 have shown that increasing viscosity of the medium induces a continuous transition from a swimming gait to a crawling gait, characterized by a decreasing undulatory frequency (Figure 7C) and an increasing curvature amplitude (Figure 7D). We asked whether our model is consistent with this load-dependent gait adaptation.

Figure 7

Download asset Open asset

Model reproduces *C.elegans* gait adaptation to external viscosity.

(A) Dark field images and the corresponding undulatory frequencies and amplitudes of adult worms (*left*) swimming in NGM buffer of viscosity 1 mPa·s, (*right*) crawling on agar gel surface. The worm head is to the right in both images. (B) Phase portrait graphs measured from worm forward movements in fluids of viscosity 10 mPa·s (blue, 3528 cycles using 50 worms), 120 mPa·s (red, 5050 cycles using 116 worms), and 5400 mPa·s (yellow, 1364 cycles using 70 worms). (**C,D**) The model predicts the dependence of undulatory frequency (C) and curvature amplitude (D) on external viscosity (black) that closely fit the corresponding experimental observations (red). (E) Phase portrait graphs predicted from the model in three different viscosities (same values as in B). Gray dashed lines indicate threshold lines for dorsoventral bending. The intersections (red circles 1, 2, 3) between the threshold line and phase portrait graphs are switch points for undulations in low, medium, high viscosity, respectively. (F) Theoretically predicted PRCs in fluids of the three different viscosities show that PRC will be shifted to the right as the viscosity of environment increases. (G) PRCs measured from optogenetic inhibition experiments in the three viscosities. Experimental PRCs were obtained using a moving average along the x-axis with $0.16 π$ in bin width and filled areas are 95% confidence interval. The tendency of shift observed in experimental PRCs verified the model prediction.

We incorporated the effect of external viscosity into our model through the bending relaxation time constant $τ_{u}$ (see Appendix). We ran our model to determine the dependence of model output on viscosity with varying viscosity $η$ . We found that model results for frequency and amplitude dependence on viscosity of the external medium are in quantitative agreement with previous experimental results (Fang-Yen et al., 2010; Figure 7C,D).

We sought to develop an intuitive understanding of how the model output changes with increasing viscosity. We recall that the model generates a proprioceptive feedback variable in the form $P = K + b \dot{K}$ (Figure 4A), and that the active muscle moment in our model undergoes a change of sign upon the proprioceptive feedback reaching either of two thresholds, $P_{t h}$ and $- P_{t h}$ . As the viscosity increases, one expects that a worm will perform a slower undulation due to the increase in external load. That is, the term $b \dot{K}$ becomes smaller. To compensate for this effect, the worm needs to undulate with a larger curvature amplitude to maintain the same level of proprioceptive feedback.

Next, we asked how the PRC depends on external viscosity. Model simulations with three different viscosities produced PRCs with similar sawtooth shape but with sharp transitions delayed in phase as the external viscosity increases (Figure 7F). We also measured PRCs from optogenetic inhibition experiments in solutions of three different viscosities (Figure 7G). Comparing the relative locations of the transitions in PRCs between the model and the data, our prediction also quantitatively agrees with the experimental results.

These results further support the model’s description of how undulatory dynamics are modulated by the external environment.

Evaluation of alternative oscillator models

Although our computational model agrees well with our experimental results, we asked whether other models could also explain our findings. We examined three alternative models based on well-known mathematical descriptions of oscillators (van der Pol, Rayleigh, and Stuart-Landau oscillators) and compared them with our original threshold-switch model and with our experimental data.

First, we tested the van der Pol oscillator, the first relaxation oscillator model (Van der Pol, 1926) which has long been applied in modeling neuronal dynamics (Fitzhugh, 1961; Nagumo et al., 1962). It is based on a second-order differential equation for a harmonic oscillator with a nonlinear, displacement-dependent damping term (see Appendix). By choosing a set of appropriate parameters, we found that the free-running waveform and phase plot of the van der Pol oscillator are highly asymmetric, but in an inverted manner (Figure 5—figure supplement 1B,F), compared with the experimental observations (Figure 2C,D). Transiently perturbing the system with the bell-shaped modulatory function over all phases within a cycle produced a similar sawtooth-shaped PRC as that observed experimentally (Figure 5—figure supplement 1N). However, the perturbed system was found to recover toward its limit cycle with a much slower rate than that of the experiments (Figure 5—figure supplement 1J). Simulations of single-side muscle inhibitions to the system produced single-sawtooth-shaped PRCs similar to those found experimentally (Figure 6—figure supplement 2B,F).

Next, we examined the Rayleigh oscillator, another relaxation oscillator model which was originally proposed to describe self-sustained acoustic vibrations such as vibrating clarinet reeds (Rayleigh, 1896). It is based on a second-order differential equation with a nonlinear, velocity-dependent damping term and it can be obtained from the van der Pol oscillator via a variable differentiation and substitution (see Appendix). From its free-running dynamics, we observed that the system exhibits a highly asymmetric waveform and phase plot that are similar to the experimental observations (Figure 5—figure supplement 1C,G). Additionally, the Rayleigh oscillator also produces similar sawtooth-shaped PRCs with respect to transient muscle inhibitions of both sides (Figure 5—figure supplement 1O), dorsal side (Figure 6—figure supplement 2C), and ventral side (Figure 6—figure supplement 2G), respectively, and system’s recovery rate after the perturbation was shown to be similar to that of the experiments (Figure 5—figure supplement 1K).

Finally, we considered the Stuart-Landau oscillator, a commonly used model for the analysis of neuronal synchrony (Acebrón et al., 2005). Its nonlinearity is based on a negative damping term which depends on the magnitude of the state variable defined in a complex domain (see Appendix). The negative damping of the system constantly neutralize the positive damping on a limit cycle, making its free-running dynamics a harmonic oscillation which shows a sinusoidal waveform (Figure 5—figure supplement 1D,H). Moreover, PRCs with respect to transient muscle inhibitions are constant with respect to phase (Figure 5—figure supplement 1P), contrary to the experiments.

We compared the results of our models with the experimental results. In the van der Pol oscillator, the free-running waveform displays a different asymmetry (Figure 5—figure supplement 1B,F) compared with the experimental observations and the perturbed system was shown to recover toward its limit cycle with a much slower rate than that of the experiments (Figure 5—figure supplement 1J). The Rayleigh oscillator reproduces a free-running waveform similar to experimental ones (Figure 5—figure supplement 1C,G) and its recovery rate toward limit cycle upon perturbation was close to that of the experiments (Figure 5—figure supplement 1K). However, its PRC (Figure 5—figure supplement 1O) showed weaker agreement with the experimental PRC compared with the threshold-switch model (Figure 5—figure supplement 1M) or the van der Pol model (Figure 5—figure supplement 1N). Of all the models tested, the threshold-switch model showed the least mean-square error with the PRC data (Figure 5—figure supplement 1M–P). We conclude that of these models, our threshold-switch model produced the best overall agreement with experiments.

We also found that two important experimental findings, the nonsinusoidal free-moving dynamics and the sawtooth-shaped PRCs can be achieved in our original model, the van der Pol and Rayleigh oscillators, which are all relaxation oscillators, but not in the Stuart-Landau oscillator, which is not a relaxation oscillator. Taken together, these results are consistent with the idea that a relaxation oscillation mechanism may underlie C. elegans motor rhythm generation.

Discussion

In this study, we used a combination of experimental and modeling approaches to probe the mechanisms underlying the C. elegans motor rhythm generation.

Our model can be compared to those previously described for C. elegans locomotion. An early model (Niebur and Erdös, 1991) assumes that a CPG located in the head initiates dorsoventral bends and that a combination of neuronal and sensory feedback mechanisms propagates the waves in the posteriorward direction. In this model, sensory feedback plays a modulatory role in producing smoother curvature waves but is not explicitly required for rhythm generation itself. Other computational models have aimed to describe how the motor circuit generates rhythmicity. Several neural models for the forward-moving circuit (Karbowski et al., 2008; Olivares et al., 2021) incorporating of all major neural components and connectivity have been developed. These models included a CPG in the head based on effective cross-inhibition between ventral and dorsal groups of interneurons. In contrast, Bryden and Cohen, 2008 developed a neural model in which each segment along the body is capable of generating oscillations. In this model, a circuit of AVB interneurons and B-type motor neurons suffices to generate robust locomotory rhythms without cross-inhibition.

Other models have examined how C. elegans adapts its undulatory wavelength, frequency, and amplitude as a gait adaptation to external load (Boyle et al., 2012; Denham et al., 2018; Izquierdo and Beer, 2018; Johnson et al., 2021). To account for these changes, these models combined the motor circuit model with additional assumptions of stretch sensitivity in motor neurons, and worm body biomechanical constraints, to create a model that reproduced the changes in undulatory wave patterns under a range of external conditions.

Previous detailed models of C. elegans locomotion have employed a relatively large number of free parameters (up to 40; Boyle et al., 2012; Karbowski et al., 2008). In our work, we sought to develop a compact phenomenological model to describe an overall mechanism of rhythm generation but not the detailed dynamics of specific circuit elements. To improve predictive power, we aimed to minimize the number of free parameters used in the model. Our model has only five free parameters, yet accurately describes a wide range of experimental findings including the nonsinusoidal dynamics of free locomotion, phase response curves to transient paralysis, and dependence of frequency and amplitude on external viscosity.

Our phase portrait analysis of worm’s free locomotory dynamics has described a previously undescribed methods for measuring the bending relaxation time scale $τ_{u}$ and the muscle moment transition time scale $τ_{m}$ (see Appendix for details), which may be compared with previous studies of worm biomechanics (Fang-Yen et al., 2010; Berri et al., 2009) and neurophysiology (Milligan et al., 1997). Fang-Yen et al., 2010 measured out a linear relationship between the bending relaxation time scale and the external viscosity by deforming the worm body in Newtonian fluids with varied viscosities in the range 1–25 mPa·s. Through an extrapolation based on that linear relationship, the relaxation time scale in 17% dextran NGM fluid (approximately 120 mPa·s in viscosity) is estimated to be $\approx 282 m s$ , which is quite close to our measured result, $τ_{u} \approx 260 m s$ . Furthermore, our measurement of the muscle moment transition time scale ( $τ_{m} \approx 100 m s$ ) is consistent with a previously measured value for muscle time scale (Milligan et al., 1997) that has been widely adopted for other detailed models of nematode locomotion (Boyle et al., 2012; Bryden and Cohen, 2008; Butler et al., 2015; Chen et al., 2011; Denham et al., 2018; Izquierdo and Beer, 2018; Johnson et al., 2021; Karbowski et al., 2008; Olivares et al., 2021; Wen et al., 2012).

In our model, the mechanism for generating rhythmic patterns can be characterized by a ‘relaxation oscillation’ process which contains two alternating sub-processes on different time scales: a long relaxation process during which the motor system varies toward an intended state due to its biomechanics under a constant active muscle moment, alternating with a rapid period during which the active muscle moment switches to an opposite state due to a proprioceptive thresholding mechanism.

The term ‘relaxation oscillation’, as first employed by van der Pol, describes a general form of self-sustained oscillatory system with intrinsic periodic relaxation/decay features (Van der Pol, 1926). The Fitzhugh-Nagumo model (Fitzhugh, 1961; Nagumo et al., 1962), a prototypical model of excitable neural systems, was originally derived by modifying the van der Pol relaxation oscillator equations. These and similar relaxation oscillators have been characterized in various dynamical systems in biology and neuroscience (Izhikevich, 2007). For example, the dynamics exhibited from the action potentials of barnacle muscles in their oscillatory modes were found to yield ‘push-pull’ relaxation oscillation characteristics (Morris and Lecar, 1981). The beating human heart was found to behave as a relaxation oscillator (Der pol b, 1940). Several studies of walking behavior in stick insects (Bässler, 1977; Cruse, 1976; Graham, 1985; Wendler, 1968) proposed that the control system for rhythmic step movements constitutes a relaxation oscillator in which the transitions between leg movements is determined by proprioceptive thresholds.

Key properties shared by these relaxation oscillators are that their oscillations greatly differ from sinusoidal oscillations and that they all consist of a certain feedback loop with a ‘discharging property’. They contain a switch component that charges an integrating component until it reaches a threshold, then discharges it again (Nave, 2007), then repeats. Many relaxation oscillators, including the van der Pol and Rayleigh models, exhibit sawtooth-shaped phase response curves (Der pol b, 1940; also see Figure 5—figure supplement 1). As shown in our experimental and model results, all the above properties have been revealed in the dynamics of C. elegans locomotive behavior, consistent with the idea that the worm’s rhythmic locomotion also results from a type of relaxation oscillator.

In our computational model, a proprioceptive component sensing the organism’s changes in posture is required to generate adaptive locomotory rhythms. What elements in the motor system could be providing this feedback? Previous studies have suggested that head and body motor neurons, including the SMDD head motor neurons and the B-type motor neurons, have proprioceptive capabilities (Wen et al., 2012; Yeon et al., 2018) and may also be involved in locomotory rhythm generation (Fouad et al., 2018; Gao et al., 2018; Kaplan et al., 2020; Xu et al., 2018). This possibility is consistent with earlier hypothesis that the long undifferentiated processes of these cholinergic neurons may function as proprioceptive sensors (White et al., 1986). In particular, recent findings (Yeon et al., 2018) have revealed that SMDD neurons directly sense head muscle stretch and regulate muscle contractions during oscillatory head bending movements.

In our model, the proprioceptive feedback variable depends on both the curvature and the rate of change of curvature. Many mechanoreceptors are sensitive primarily to time derivatives of mechanical strain rather than strain itself; for example, the C. elegans touch receptor cells exhibit such a dependence (Eastwood et al., 2015; O'Hagan et al., 2005). The ability of mechanosensors to sense the rate of change in C. elegans curvature has been proposed in an earlier study (Butler et al., 2015) in which it was hypothesized that the B-type motor neurons might function as a proprioceptive component in this manner. Mechanosensors encoding a simultaneous combination of deformation and velocity have been observed in mammalian systems including rapidly-adapting (RA) and intermediate-adapting (IA) sensors in the rat dorsal root ganglia (Rugiero et al., 2010). Proprioceptive feedback that involves a linear combination of muscle length and velocity was also suggested by a study of C. elegans muscle dynamics during swimming, crawling, and intermediate forms of locomotion (Butler et al., 2015). In our phenomenological model, the motor neuron constituent may represent a collection of neurons involved in motor rhythm generation. Therefore, the proprioceptive function posited by our model might also arise as a collective behavior of curvature-sensing and curvature-rate-sensing neurons.

Further identification of the neuronal substrates for proprioceptive feedback may be possible through physiological studies of neuron and muscle activity using calcium or voltage indicators. Studies of the effect of targeted lesions and genetic mutations on the phase response curves will also help elucidate roles of specific neuromuscular components within locomotor rhythm generation.

In summary, our work describes the dynamics of the C. elegans locomotor system as a relaxation oscillation mechanism. Our model of rhythm generation mechanism followed from a quantitative characterization of free behavior and response to external disturbance, information closely linked to the structure of the animal’s motor system (Gutkin et al., 2005; Nadim et al., 2012; Schultheiss et al., 2011; Smeal et al., 2010). Our findings represent an important step toward an integrative understanding of how neural and muscle activity, sensory feedback control, and biomechanical constraints generate locomotion.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Strain, strain background (E. coli)	OP50	CGC	Fang-Yen Lab Strain Collection: OP50 RRID:WB-STRAIN:WBStrain00041971	OP50
Strain, strain background (C. elegans)	YX148	Fouad et al., 2018	Fang-Yen Lab Strain Collection: YX148	qhIs1[Pmyo-3::NpHR::eCFP; lin-15(+)]; qhIs4[Pacr-2::wCherry]
Strain, strain background (C. elegans)	YX119	Fouad et al., 2018	Fang-Yen Lab Strain Collection: YX119	qhIs1[Pmyo-3::NpHR::eCFP; lin-15(+)]; unc-49(e407)
Strain, strain background (C. elegans)	YX205	Leifer et al., 2011	Fang-Yen Lab Strain Collection: YX205	hpIs178[Punc-17::NpHR::eCFP; lin-15(+)]
Strain, strain background (C. elegans)	WEN001	Fouad et al., 2018	Fang-Yen Lab Strain Collection: WEN001	wenIs001[Pacr-5::Arch::mCherry; lin-15(+)]

Worm strains and cultivation

Request a detailed protocol

C. elegans were cultivated on NGM plates with Escherichia coli strain OP50 at 20°C using standard methods (Sulston and Hodgkin, 1988). Strains used and the procedures for optogenetic experiments are described in the Key resources table and Appendix. Preparation of OP50 and OP50-ATR plates were as previously described (Fouad et al., 2018). All experiments were performed with young adult (< 1 day) hermaphrodites synchronized by hypochlorite bleaching.

Locomotion and phase response analyses

Request a detailed protocol

To perform quantitative recordings of worm behavior, we used a custom-built optogenetic targeting system as previously described (Fouad et al., 2018; Leifer et al., 2011). Analysis of images for worm’s body posture was performed using a previously developed custom software (Fouad et al., 2018). The anterior curvature is defined as the average of the curvature over body coordinate 0.1–0.3; excluding the range from 0 to 0.1 avoided measurement of high-frequency movements of the worm’s anterior tip. Descriptions of the apparatus and image analyses are available in Appendix.

For phase response experiments, opsin-expressing worms were illuminated using a brief laser pulse (532 nm wavelength, 0.1 or 0.055 s duration, irradiance 16 mW/mm²) in the head region (0–0.25 body coordinate). A total of 10 trials with 6 s interval between successive pulses were performed for each animal. Trials in which the worms did not maintain forward locomotion were censored. To generate the phase response curve (PRC), we calculated the phase of inhibition of each trial and the resulting phase shift. Details of calculations for the averaged PRC are given in Appendix.

All the data and image analysis codes used in the manuscript are available at Dryad (archived at https://doi.org/10.5061/dryad.wwpzgmsk2).

Computational modeling

Request a detailed protocol

Our primary model is based on a novel neural control mechanism incorporated with a previously described biomechanical framework (Fang-Yen et al., 2010; Gray and Lissmann, 1964; Wallace, 1968). A proprioceptive signal is defined by a linear combination of bending curvature and rate of change of curvature. When the signal reaches a threshold, a switching command is initiated to reverse the direction of muscle moment. The worm’s curvature relaxes toward the opposite direction, and the process repeats, creating a dorsoventral alternation. Detailed descriptions of implementation and fitting procedure of this model and alternative models are available in Appendix. All codes for modeling analyses are available at Dryad (https://doi.org/10.5061/dryad.wwpzgmsk2).

Appendix 1

Strains and plates preparations

Transgenic strains used in this study are listed as follows: YX148 qhIs1[Pmyo-3::NpHR::eCFP;lin-15(+)];qhIs4[Pacr-2::wCherry], YX119 unc-49(e407);qhIs1[Pmyo-3::NpHR::eCFP;lin-15(+)], YX205 hpIs178[Punc-17::NpHR::eCFP;lin-15(+)], WEN001 wenIs001[Pacr-5::Arch::mCherry;lin-15(+)].

For optogenetic experiments, worms were cultivated in darkness on plates with OP50 containing the cofactor all-trans retinal (ATR). For control experiments and free-moving experiments, worms were cultivated on regular OP50 NGM plates without ATR. To make OP50-ATR plates, we added 2 µL of a 100 mM solution of ATR in ethanol to an overnight culture of 250 µL OP50 in LB medium and used this mixture to seed 6 cm NGM plates.

Locomotion analysis

To analyze worm locomotion in viscous fluids, we placed worms in dextran solutions in chambers formed by a glass slide and a coverslip separated by 125-µm-thick polyester shims (McMaster-Carr 9513K42). For viscosity-dependence experiments, we used 5%, 17%, and 35% (by mass) solutions of dextran (Sigma-Aldrich D5376, average molecular weight 1,500,000–2,800,000) in NGMB. These solutions were measured to have viscosities of 10, 120, and 5400 mPa·s (Fang-Yen et al., 2010), respectively. We used a 17% dextran solution for all other experiments. NGMB consists of the same components as NGM media (Stiernagle, 2006), but without agar, peptone, or cholesterol.

We recorded image sequences using a custom-built optogenetic targeting system based on a Leica DMI4000B microscope under 10X magnification with dark field illumination provided by red LEDs. Worm images were recorded at 40 Hz with an sCMOS camera (Photometrics optiMOS). We used a custom-written C++ software (Fouad et al., 2018) to perform real-time segmentation of the worm during image acquisition. The worm was identified in each image by its boundary and centerline calculated from a binary image. Anterior-posterior orientation was noted visually during the recording. Segmentation information, including coordinates of the worm boundary and centerline, was saved to disk along with the corresponding image sequences.

Post-acquisition image analysis was performed using a custom MATLAB (Mathworks) similar to previous reports (Fouad et al., 2018). The worm centerline of each image was smoothed using a cubic spline fit. We calculated curvature $κ$ as the dot product between the unit normal vector to the centerline and derivative of the unit tangent vector to the centerline with respect to the body coordinate. Dimensionless curvature $K$ was calculated as the product of $κ$ and the worm body length $L$ represented by length of the centerline. Since the segmentation was relatively noisy at the tips of the worm, we excluded curvature in the anterior and posterior 5% of the body length. The worm’s direction of motion was identified by calculating the gradients in the curvature over time and body coordinate, and image sequences in which the worm performed consistent forward movement (lasting at least 4 s) were selected for analysis. The anterior curvature $K (t)$ was defined as the average of the dimensionless curvature over body coordinate 0.1-0.3; this range avoided high-frequency movements of the anterior tip of the animal.

To quantify oscillatory dynamics during forward locomotion, we identified undulatory cycles from the time sequence of anterior curvature in each worm. Local extrema along each sequence were identified and portions between consecutive local maxima were defined as individual cycles. To minimize the effects of changes in the worm’s frequency, we excluded cycles whose period deviated by more than 20% from the average period of all worms’ undulations in each experimental session.

For the ease of computing average dynamics, we converted individual cycles from a time-dependent to a phase-dependent curvature by uniformly rescaling each cycle to a phase range of $2 π$ . The averaged curvature within one cycle was then computed by averaging all individual cycles in the phase domain: $⟨ K (ϕ) ⟩ = \sum_{i = 1}^{N} K_{i} (ϕ) / N$ . Similarly, the averaged phase derivative of curvature within one cycle was calculated as $⟨d K / d ϕ⟩ = \sum_{i = 1}^{N} (d K_{i} / d ϕ) / N$ .

Stability of the worm’s head oscillation

To examine the stability of the worm’s head oscillation during forward locomotion, we analyzed head oscillations of worms that were optogenetically perturbed with 0.1 s muscle inhibitions and estimated their recovery dynamics after being deviated from the normal oscillation due to the perturbation.

To illustrate the oscillation dynamics, we use a two-dimensional variable, $x = (K, ξ \dot{K})$ in the unit of curvature where $ξ = 0.135 s$ is a scaling factor. In Figure 3—figure supplement 2, we depicted the closed trajectory (black) in the plane spanned by the variables $K$ and $ξ \dot{K}$ for the head oscillation of unperturbed moving worms (this coordinate plane is in fact a linearly scaled version of the phase plane spanned by the variables $K$ and $\dot{K}$ ), which we call as the normal cycle of the worm’s head oscillation.

Next, we defined an amplitude variable $d$ that represents the normalized deviation to the normal cycle. If the oscillator is stable, the closed orbit for the unperturbed dynamics is usually called the stable limit cycle. Here, we stick to the notion of normal cycle instead of using ‘limit cycle’ to avoid any confusion on the stability of the worm’s head oscillation. For any phase state of an individual oscillation, the normalized deviation to the normal cycle is defined as $d (ϕ) = (D (ϕ) - D^{C} (ϕ)) / D^{C} (ϕ)$ . Here, $D (ϕ)$ is distance to the center of oscillation on the phase plane, which is set to the origin, such that $D (ϕ) = \sqrt{K {(ϕ)}^{2} + {(ξ \dot{K} (ϕ))}^{2}}$ where $ϕ$ denotes the phase value of the current state estimated by the four-quadrant inverse tangent of the variable pair $(K, ξ \dot{K})$ . In this expression, $D^{C} (ϕ)$ denotes the distance to the center of oscillation that is evaluated exactly on the normal cycle at phase $ϕ$ .

Using the deviation to the normal cycle to describe the amplitude of the worm’s head oscillation, we collected the amplitude dynamics over time for all periods of the worm’s head oscillations during which no illumination pulse occurs, that is, all periods of locomotion between two consecutive illumination pulses. We grouped the amplitude dynamics into bins according to their initial amplitudes and then calculated the collective amplitude dynamics for each bin. As shown in Figure 3—figure supplement 1, the collective amplitude variable $d$ converges to zero after roughly 0.5 s regardless of the initial amplitude. This result indicates that the worm’s head oscillation returns to its normal oscillation after being perturbed and that the normal cycle may represent a stable limit cycle for the oscillation.

Phase isochron map and vector field for the worm’s head oscillation

On the normal cycle we define the phase of the oscillation as $ϕ^{C} (t) = ω_{0} \cdot t_{m o d T_{0}}$ , where $ω_{0} = 2 π / T_{0}$ is the angular frequency of normal oscillation (the calculation of $T_{0}$ will be described in the next subsection) and we determined the initial phase ( $ϕ^{C} = 0$ ) to be when $K$ reaches local maximum (or $x = (K_{m a x}, 0)$ ) and hence $ϕ^{C} = π$ to be when $K$ reaches local minimum (or $x = (K_{m i n}, 0)$ ). In this way, we parameterized the normal cycle by defining a bijective map between phases and state points $Φ (x^{C}) = ϕ^{C}$ .

The map $Φ (x) = ϕ$ has been well defined for all the state points on the normal cycle $C$ . We next estimate the phases for points off the normal cycle. By definition (Izhikevich, 2007), if $x_{0}$ is a point on the normal cycle and $y_{0}$ is a point off the normal cycle, then $y_{0}$ will have the same phase as $x_{0}$ if the trajectory starting at the initial point $y_{0}$ off the normal cycle converges to the trajectory starting at the initial point $x_{0}$ on the normal cycle as time goes to infinity. Here, we define the set of all state points off the normal cycle having the same phase as a point $x_{0}$ on the normal cycle as the isochron (Winfree, 2001) for phase $ϕ_{0} = Φ (x_{0})$ .

In our analysis, it was not possible to define an isochron according to the theoretical definition since data were always recorded in a finite time period during experiments. We used an alternative way to estimate all isochrons on the phase plane for the worm’s head oscillation. For each individual trial of illumination, we observed that, due to the optogenetic inhibition, the variable $\dot{K}$ quickly decayed toward zero value immediately after the illumination and then recovered after approximately 0.3 s as the oscillation converged to a normal oscillation. Therefore, by finding the local minimal of $\dot{K}$ immediately after each illumination pulse, we located the point at which the paralyzing effect is just removed and after which the oscillation starts a free resumption to normal oscillation. We call this point the 'notch point' $x_{N}$ as it can be clearly seen from the phase plot (as shown in Figure 3E, H and K). After the notch point $x_{N}$ , the oscillation will proceed to its next phase states $x (ϕ = 0)$ and $x (ϕ = π)$ (or vice-versa), both of which can be easily identified through peak finding from the curvature dynamics $K$ . Hence, we obtained two sub-trajectories from the oscillation: one $x_{N} \to x (ϕ = 0)$ , and the other $x_{N} \to x (ϕ = π)$ . Next to determining the timing of the notch point $t (x_{N})$ , we determined the phase of the notch point in the following steps: (1) we computed the phase of the state on the normal cycle at time $t (x_{N})$ as if the perturbation had not occurred, which is $ϕ_{u} (x_{N}^{C}) = ω_{0} {(t (x_{N}) - t (x (ϕ = 0)))}_{m o d T_{0}}$ , or $ϕ_{l} (x_{N}^{C}) = ω_{0} {(t (x_{N}) - t (x (ϕ = π)))}_{m o d T_{0}} - π$ . Here, $ϕ (x_{N}^{C})$ was computed twice using phase states $x (ϕ = 0)$ [subscripted with $u$ ] and $x (ϕ = π)$ [subscripted with $l$ ] as references, respectively; (2) we calculated the induced phase shift $P R C (t_{i l l u m})$ and the phase of the notch point is $ϕ (x_{N}) = ϕ (x_{N}^{C}) - P R C (t_{i l l u m})$ . After obtaining the sub-trajectories $x_{N} \to x (ϕ = 0)$ and $x_{N} \to x (ϕ = π)$ and calculating the phase of $x_{N}$ , we then estimated the phase values for all the points within each of the two sub-trajectories through linear interpolation.

Following the above steps, we calculated the phase values for all the state points on the phase plane that have been recorded from the optogenetic experiments. We then applied a 2-D moving average (using the angular statistics method) for the obtained phase values over the phase plane to smooth out the isochron map. Finally, we used a linear 2-D interpolation to obtain a phase isochron map with a finer resolution as shown in Figure 3—figure supplement 2.

To compute the vector field of the worm’s head oscillation, we collected all the sub-trajectories $x_{N} \to x (ϕ = 0)$ and $x_{N} \to x (ϕ = π)$ that are defined above and took derivative of the trajectories with respect to time. Thus, by collecting all the phase states $(K, ξ \dot{K})$ and their corresponding time derivatives $(d K / d t, d (ξ \dot{K}) / d t)$ that describe the tangent vectors of trajectories, we generated the raw form of vector field for the worm’s head oscillation. Again, we applied a 2-D moving average for the raw outcome over the phase plane to smooth out the vector field. We used a linear 2-D interpolation to obtain a vector field with an appropriate number of quivers to be displayed (Figure 3—figure supplement 2).

Phase response analysis

To generate phase response curves (PRCs) from optogenetic inhibition experiments, each trial’s illumination phase $ϕ$ , as well as the induced phase shift $F$ , were calculated. To calculate the two variables, the animal’s phase of oscillation was estimated based on timings of local extrema identified from the time-varying curvature profiles via a peak finding method. Specifically, (i) the occurrence of illumination of the trial was set to $t = T_{i l l u m}$ ; $t = 0$ was set at the beginning of each experiment. (ii) Around the illumination, timings of the two local maxima of curvature immediately before and after were identified as the two zero-phase points of the oscillation before and after the illumination, respectively. Here, the timings are denoted as $T Z_{- 2}$ , $T Z_{- 1}$ , $T Z_{+ 1}$ , and $T Z_{+ 2}$ , in the ascending order of time. (iii) Similarly, timings of the two local minima of curvature immediately before and after the illumination were identified as the two half-cycle-phase points before and after the illumination, respectively. Here, the timings are denoted as $T H_{- 2}$ , $T H_{- 1}$ , $T H_{+ 1}$ , and $T H_{+ 2}$ , in the ascending order of time. (iv) With these measurements, cycle period $T_{0}$ was computed as $T_{0} = (T Z_{+ 2} - T Z_{+ 1} + T Z_{- 1} - T Z_{- 2} + T H_{+ 2} - T H_{+ 1} + T H_{- 1} - T H_{- 2}) / 4$ , so the angular frequency of undulation $ω_{0} = 2 π / T_{0}$ ( $T_{0}$ was computed as the average of differences of adjacent local maxima/minima before and after illumination; multiple cycles were used here to reduce noise). In addition, the illumination phase $ϕ$ of each individual trial was computed as $ϕ_{u} = ω_{0} {(T_{i l l u m} - T Z_{- 1})}_{m o d T_{0}}$ , $ϕ_{l} = ω_{0} {(T_{i l l u m} - T H_{- 1} + T_{0} / 2)}_{m o d T_{0}}$ , and the corresponding phase shift $F$ was computed as $F_{u} = ω_{0} {(T Z_{+ 1} - T Z_{- 1})}_{m o d T_{0}} - π$ , $F_{l} = ω_{0} {(T H_{+ 1} - T H_{- 1} + T_{0} / 2)}_{m o d T_{0}} - π$ . Here, phase of illumination and the corresponding phase shift were computed twice using zero [subscripted with $u$ ] and half-cycle [subscripted with $l$ ] phase points as references, respectively.

We generated 2-D scatter plots for all trials with illumination phase as x coordinate and the corresponding phase shift as y coordinate. To visualize the distribution of the scatter points we generated bivariate histogram plots by grouping the data points into 2-D bins with 25 bins on both dimensions covering the range $[0, 2 π]$ for x and range $[- π, π]$ for y. To indicate average tendency of phase shift depending on phase of illumination, we calculated a mean-curve representation of PRCs via a moving average operation. In this process, each mean was calculated over a sliding window of width $0.16 π$ along the direction of $ϕ$ from 0 to $2 π$ . The 95% confidence interval relative to each window of data points was also computed and an integral number of them were displayed as filled area around the PRC. Through the computation, all statistical calculations followed the rules of directional statistics (Fisher et al., 1993) since $ϕ$ and $F$ are circular variables defined in radians.

Phase response curves from perturbations of other body regions

We asked how phase responses for the other regions of the body would compare to that of the anterior region. We conducted optogenetic experiments that inhibited Pmyo-3::NpHR transgenic worms by transiently illuminating 0.1–0.3 (anterior), 0.4–0.6 (middle), and 0.6–0.8 (posterior) of the body length, respectively. We found that the amplitude of the sawtooth feature of PRC tends to decrease as the perturbation occurs further from the head (Figure 3—figure supplement 5A,E,I). We also noticed that, for the same perturbed region, the PRC shape remains unaffected regardless of the region at which the dynamics were analyzed (see Figure 3—figure supplement 5A–C,D–F,G–I, respectively), suggesting that posterior regions of a freely moving worm follow their anterior neighbors with a constant phase offset. Taken together, these results suggest that a main rhythm generator may operate near the head of the worm to produce primary oscillations during forward locomotion. The sawtooth-shape feature of the PRC becomes stronger if the perturbation hits closer to the rhythm generator (Figure 3—figure supplement 5A) or becomes weaker if the perturbation indirectly affects it (Figure 3—figure supplement 5E,I)

The relaxation oscillator model for locomotor wave generation

We first developed a relaxation oscillator model to simulate the rhythm generation during C. elegans forward locomotion. In this model, we incorporated a novel neuromuscular mechanism with a previously described biomechanical framework (Fang-Yen et al., 2010). Here, we only simulated the bending rhythms generated from the head region; the wave propagation dynamic is out of the scope of our study. Our phenomenological model does not contain detailed activities of individual neurons but focus on describing key neuromuscular mechanisms and their contributions to the rhythm generation.

To produce model variables that can be directly compared with experimental observations of moving animals, a biomechanical framework was first developed to describe worm’s behavioral dynamics in its external environments. Following previous derivations for C. elegans biomechanics (Fang-Yen et al., 2010), the relationship between animal behavioral outputs and the active muscle moments can be described as follows:

\begin{matrix} C_{N} \frac{\partial y}{\partial t} + a \frac{\partial^{2} κ}{\partial s^{2}} + a_{v} \frac{\partial}{\partial t} (\frac{\partial^{2} κ}{\partial s^{2}}) = m_{a} . \end{matrix}

In Equation A1, the first term indicates the external viscous force that is transverse to the body segment where $C_{N}$ is the coefficient of viscous drag to the transverse movement and $y$ denotes the lateral displacement of body segment; the second term indicates the internal elastic force where $a$ is the bending modulus of the worm body; the third term indicates the internal viscous force where $a_{v}$ is the coefficient of the body internal viscosity. On the right side of Equation A1 is the active muscle moment $m_{a}$ .

Taking the second partial derivative with respect to body coordinate $s$ on both sides of Equation A1 and, using the linear relation under the small-amplitude approximation, $κ \approx y_{s s}$ , we arrive at:

\begin{matrix} C_{N} \frac{\partial κ}{\partial t} + a \frac{\partial^{4} κ}{\partial s^{4}} + a_{v} \frac{\partial}{\partial t} (\frac{\partial^{4} κ}{\partial s^{4}}) = \frac{\partial^{2} m_{a}}{\partial s^{2}} . \end{matrix}

Under the assumptions of small-amplitude undulations and a fixed wavelength $λ$ down the worm body, $κ$ can be considered as a travelling sinusoidal wave with a small deviation, $κ (s, t) = κ_{0} \sin (2 π s / λ - ω t) + δ$ , which leads to an approximation, $κ_{s s s s} \approx {(2 π / λ)}^{4} κ$ . Plugging this approximation into Equation A2 while keeping $s$ fixed, after some rearrangement, one gets:

\begin{matrix} κ + \frac{C_{N} {(\frac{λ}{2 π})}^{4} + a_{v}}{a} \dot{κ} = \frac{λ^{4}}{{(2 π)}^{4} a} \frac{\partial^{2} m_{a}}{\partial s^{2}} . \end{matrix}

In terms of the dimensionless curvature $K = κ \cdot L$ and dimensionless muscle moment

\begin{matrix} M_{a} = \frac{λ^{4} L}{{(2 π)}^{4} a} \frac{\partial^{2} m_{a}}{\partial s^{2}}, \end{matrix}

we can rewrite Equation A3 as:

\begin{matrix} K + τ_{u} \dot{K} = M_{a}, \end{matrix}

where

\begin{matrix} τ_{u} = \frac{C_{N} {(\frac{λ}{2 π})}^{4} + a_{v}}{a}, \end{matrix}

and we note that Equations A5 and A6 yield Equation 1. In Equation A6, both the wavelength $λ$ and the normal viscous drag coefficient $C_{N}$ vary with the fluid viscosity $η$ (Berri et al., 2009 ; Fang-Yen et al., 2010).

The above biomechanical framework in our model treats the worm’s body segment as a viscoelastic rod and describes how the body segment bends under the forces provided by the active muscle moment. However, the simulated oscillation in $K$ comes from the rhythmicity of the active muscle moment which originates from the hypothesized neuromuscular mechanism described by the following relaxation-oscillation process:

Proprioceptive feedback is sensed as a linear combination of the current curvature value and the current rate of change of curvature, $P = K + b \dot{K}$ (black curve in Figure 4D).
During movement of bending, this proprioceptive feedback is constantly compared with two threshold values $P_{t h}$ and $- P_{t h}$ (gray dashed bars in Figure 4D).
Once the feedback reaches either of the thresholds (the switch points as indicated by red circles in Figure 4D), a switch command is initiated (blue square wave in Figure 4E).
The switch command triggers the active muscle moment to change toward the opposite saturation value (black curve in Figure 4E).

To simulate the switch-triggered muscle transition we used a modified logistic function: $M_{a} (t) = \pm M_{0} ∙ \tanh (t / 2 τ_{m})$ . Here, the plus sign indicates the dorsal-to-ventral muscle moment transition while the minus sign indicates the opposite direction.

To initiate the oscillation in our model we set the system to bend toward the ventral side by setting $M_{a} |_{t = 0} = M_{0}$ and $K |_{t = 0} = 0$ . During forward locomotion, the active muscle moment oscillates by undergoing a relaxation oscillation process: a relaxation subperiod during which $M_{a}$ stays at a saturated bending state ( $M_{0}$ for ventral bending, $- M_{0}$ for dorsal bending), alternating between a shorter subperiod during which $M_{a}$ quickly transits toward the opposite state due to effects described in iii and iv. The bending curvature $K (t)$ which is driven by $M_{a}$ in an exponential decaying manner (Equation A5) follows the rhythmic activity of $M_{a}$ , thereby also exhibiting an oscillatory dynamic (Figure 4B).

This relaxation oscillator model reproduces two key features of free locomotion that we observed from experiments. First, freely moving worms exhibit nonsinusoidal curvature waveform with an intrinsic asymmetry: bending toward the ventral or dorsal directions occurs slower than straightening toward a straight posture during each locomotory cycle (Figure 4F). Second, dynamic of the active muscle moment shows a trapezoidal waveform during forward locomotion (Figure 2D Inset and Figure 4E). These results are independent of external conditions but reflect intrinsic properties of the neuromuscular mechanisms underlying locomotion rhythm generation.

Note that parameters $M_{0}$ , $τ_{u}$ , and $τ_{m}$ were estimated from data of free locomotion using phase portrait techniques described in the following subsection. Parameters $b$ and $P_{t h}$ were yet degenerate in this model of free locomotion. Here, we temporarily set $b = 0$ and then set $P_{t h}$ such that the oscillatory period predicted by model matched the average period measured from experiments with a minimum squared error:

\begin{matrix} P_{t h} = \underset{P_{t h} > 0}{a r g m i n} {(T_{m o d e l} (P_{t h}) - T_{e x p e r i m e n t})}^{2} . \end{matrix}

The nondegeneracy of $b$ and $P_{t h}$ was determined by fitting the model to the experimental PRC as described in the later subsection so that all the parameters for the model are provided as $M_{0} = 8.45$ , $τ_{u} = 260 m s$ , $τ_{m} = 100 m s$ , $b = 46 m s$ , and $P_{t h} = 2.33$ .

Measuring bending relaxation time scale and amplitude of active muscle moment

To estimate these two parameters, we applied a heuristic method that uses the shape properties of C. elegans free-running phase plot (Figure 2D). From the curve in the figure, we noticed two ‘flat’ portions symmetrically distributed at quadrant I and III on the phase plane. Recalling Equation 1 (or Equation A5): $K + τ_{u} \dot{K} = M_{a} (t)$ , the two flat regions indicate that the scaled active muscle moment, $M_{a} (t)$ , is nearly constant during the corresponding time bouts.

We then computed the linear correlation between variables $K$ and $\dot{K}$ to identify the two ‘flat’ regions and, through linear fits, obtained two linear relations respectively: $⟨ K ⟩ + τ_{1} ⟨ \dot{K} ⟩ = M_{1}$ (region 1) and $⟨ K ⟩ + τ_{2} ⟨ \dot{K} ⟩ = M_{2}$ (region 2). Thus, the bending relaxation time scale $τ_{u}$ and the amplitude of the scaled active muscle moment are estimated as $\hat{τ_{u}} = (τ_{1} + τ_{2}) / 2$ and ${\hat{M}}_{0} = (| M_{1} | + | M_{2} |) / 2$ , respectively.

The above method used the phase plot measured from locomotion of worms swimming in a 17% dextran solution (120 mPa·s viscosity) as an example. However, it is also valid for estimating parameters of locomotion in other viscosities.

Measuring active moment transition time scale

With $τ_{u}$ (estimated from the above method), $⟨K⟩$ and $⟨ \dot{K} ⟩$ (measured from locomotion) plugged to the left side of Equation 1, we were able to compute the waveform of the scaled active muscle moment $M_{a} (t)$ on the right side of Equation 1. As expected and shown in Figure 2D Inset, the curve of $M_{a} (t)$ is roughly centrally symmetric around point $(T_{0} / 2, 0)$ on the plane, with two plateau portions indicating two saturated states for dorsal and ventral muscle contractions, respectively.

Between the two plateau portions represents a period during which the active muscle moment is undergoing a ventral-to-dorsal (or vice-versa) transition. We used a modified logistic function to model the ventral-to-dorsal muscle moment transition (substituting $t$ with $- t$ for transition in the other direction):

\begin{matrix} M_{a} (t) = M_{0} \cdot \tanh (\frac{t}{τ_{m}}) . \end{matrix}

To estimate $τ_{m}$ , the exponential time constant for the transition of active muscle moment, we took the time derivative of Equation A8 and took the absolute value of the resultant:

\begin{matrix} | \frac{d M_{a}}{d t} | = \frac{M_{0}}{τ_{m}} \cdot \frac{\exp (2 t / τ_{m})}{{(1 + \exp (2 t / τ_{m}))}^{2}} . \end{matrix}

We notice that when $t = 0$ , the maximum of $| d M_{a} / d t |$ is achieved and the value is $M_{0} / τ_{m}$ . On the other hand, the maximum of $|d M_{a} / d t|$ can be obtained from the experimental observation by simply finding the peak of | $d M_{a} / d t$ | curve where $M_{a} = ⟨K (t)⟩ + \hat{τ_{u}} \cdot ⟨d K (t) / d t⟩$ . Thus, $τ_{m}$ can be estimated as:

\begin{matrix} τ_{m} = \overset{{\hat{M}}_{0}}{} \cdot {| \frac{d M_{a}}{d t} |}_{m a x}^{- 1} . \end{matrix}

Parameter estimation

For our original threshold-switch model, parameters $τ_{u}$ , $τ_{m}$ , and $M_{0}$ were estimated from free locomotion experiments as described above. These three parameters nearly fully determine the biomechanical framework of C. elegans bending movements (governed by Equation A5 and A8). On the other hand, parameters $b$ and $P_{t h}$ describe the proprioceptive feedback and the threshold-switch features in our model. Specifically, they characterize two threshold lines, $K + b \dot{K} = \pm P_{t h}$ (as shown in Figure 4C). The two switch points—defined by the intersection between the phase trajectory and the threshold lines on the phase plane—determine the timing of switches for the active muscle moment (see Figure 4C–E). We noted that the model behavioral output of free locomotion is degenerate with respect to these two parameters; the same outcome would be produced if the threshold lines cross the same pair of switch points. To first determine the free-moving dynamic as well as the switch points, we temporarily set $b = 0$ and then set $P_{t h}$ such that the oscillatory period defined by model matched the average period measured from the experiments.

To obtain the nondegeneracy of $P_{t h}$ and $b$ , we fit our model to the experimental phase response curve using a global optimization procedure. Full procedure for the determination of $b$ and $P_{t h}$ is given below.

Modeling worm oscillations in varied environments

Differences in various environments will change only those parameters that are related to contact with external forces whereas parameters related to oscillator’s internal properties will not be affected. In terms of the internal parameters of our model, we used values that were previously determined, which are $τ_{m} = 100 m s$ , $M_{0} = 8.45$ , $b = 46 m s$ , $P_{t h} = 2.33$ . For the exogenous parameters, only the time constant of undulation, $τ_{u}$ , varies according to external conditions. According to Equation A6, $τ_{u}$ is explicitly determined in terms of other physical parameters, including biomechanical parameters measured in previous work (Fang-Yen et al., 2010): the internal viscosity of worm body is measured as $a_{v} = 5 \cdot 10^{- 16} N m^{3} s$ ; the bending modulus of worm body is measured as $a = 9.5 \cdot 10^{- 14} N m^{3}$ ; $C_{N} = 31 η$ is the coefficient of viscous drag for movement normal to the body (Katz et al., 1975), where $η$ is the fluid viscosity. According to previous measurements of undulatory wavelengths in different viscous solutions (Fang-Yen et al., 2010), we applied a logarithmic fit to the data points, yielding $λ / L = - 0.158 \log_{10} (η / η_{0}) + 1.5$ for a continuous model realization in undulatory frequency and amplitude. Here, $λ$ is the wavelength and $η_{0} = 1 m P a \cdot s$ .

Alternative models for locomotor wave generation

To further evaluate the performance of our original model, we explored three alternative models for simulating locomotory rhythm generation to make comparisons across these models and the experimental observations. Alternative models are based on three previously studied self-oscillator models described by 2-D nonlinear systems: the van der Pol, Rayleigh, and Stuart-Landau oscillators.

First, we developed a model oscillator in the form taken from the van der Pol Oscillator:

\begin{matrix} \ddot{K} + a_{V} (b_{V} K^{2} - 1) \dot{K} + ω_{V}^{2} K = 0, \end{matrix}

where $K$ indicates the nondimensional bending curvature. This model has a nonlinear damping term with a coefficient depending on $K$ . Simulated oscillation is a limit cycle of the model (Figure 5—figure supplement 1B,F), given parameters $a_{V} = - 0.026 s^{- 1}, b_{V} = - 2.04, ω_{V} = 5.51 s^{- 1} .$

Second, we developed a model oscillator by modifying the Rayleigh Oscillator:

\begin{matrix} \ddot{K} + a_{R} (b_{R} {\dot{K}}^{2} - 1) \dot{K} + ω_{R}^{2} K = 0, \end{matrix}

where $K$ again indicates the nondimensional bending curvature. This model has a nonlinear damping term with a coefficient depending on $\dot{K}$ . Simulated oscillation is a limit cycle of the model (Figure 5—figure supplement 1C,G), given parameters $a_{R} = 2.73 s^{- 1}, b_{R} = 0.0023 s^{2}, ω_{R} = 5.6 s^{- 1} .$

Third, we developed a model oscillator by modifying the Stuart-Landau Oscillator:

\begin{matrix} \dot{Z} + (\frac{l}{2} {| Z |}^{2} - σ) Z = 0. \end{matrix}

Here, the system is described in the complex domain where $Z = Z_{r} + i Z_{i}$ , $l = l_{r} + i l_{i}$ are complex variables, and $σ$ is real. We let $Z_{r}$ , the real part of $Z$ , denote the nondimensional curvature $K$ . This model has a nonlinear damping term with coefficient depending on $| Z |$ . Simulated oscillation is a limit cycle of the model (Figure 5—figure supplement 1D,H), given parameters $l_{r} = 0.54 s^{- 1}, l_{i} = 0.52 s^{- 1}, σ = 5.54 s^{- 1} .$

The three alternative models produce free-running oscillatory dynamics with similar frequency and amplitude as measured from worms swimming in fluids with viscosity $120 m P a \cdot s$ .

Simulation of optogenetic inhibition

According to our experimental observations on the effect of the optogenetic muscle inhibition (Figure 3A,B), paralysis of muscles of the illuminated region initiated upon illumination (defined as $t = 0$ for Figure 3B) and reached maximal effect approximately at $t = 0.3 s$ . Here, we modeled the process of muscle inhibition by multiplying the scaled active muscle moment, $M_{a}$ , with a factor, $1 - Q (Δ t)$ , as a function of the time interval $Δ t$ in a bell-shaped form (Figure 4—figure supplement 1, Equation A14).

As described in our model, the dorsoventrally alternating feature of the active muscle moment during locomotion are described by the dynamics of $M_{a} (t)$ . Specifically, $M_{a} (t)$ is positive when ventral muscles contract and dorsal muscles relax, and negative for the other half of the cycle. Therefore, in our threshold-switch model, specifically inhibiting dorsal- or ventral- or both-side muscles was computationally equivalent to conditionally modulating $M_{a} (t)$ with the bell-shaped modulating function depending on the sign of $M_{a} (t)$ .

For simulating inhibition process in the three alternative models, we factored out a specific term from individual model equations as a generalized active muscle moment. We applied the bell-shaped modulating function to this term conditionally for each individual model. Detailed descriptions of implementing modeled inhibitions in alternative models are available from below.

To get a deeper understanding of how phase response curves are related to systems dynamics during wave generation, we systematically simulated transient muscle inhibitions on individual model oscillators at different times within a cycle period to generate model PRCs. To do that, we theoretically simulated the process of muscle inhibition by multiplying model active muscle moment with a modulatory factor, $1 - Q (Δ t)$ , which has a bell-shaped profile (Figure 4—figure supplement 1):

\begin{matrix} Q (Δ t) = \frac{H}{(1 + {| \frac{Δ t - r}{p} |}^{2 q})}, \end{matrix}

where $r = 0.3 s$ is the timing of the occurrence of maximal paralysis according to our experimental observations on the effect of muscle inhibition (Figure 3A,B), $H$ indicated the maximal degree of paralysis, and $p$ , $q$ measure the paralyzing rate and duration, respectively. To ensure sufficient smoothness during computation, we let $p = 0.3 \cdot 10^{- 1 / q}$ so that $Q |_{Δ t = 0} > 0.99$ . Note that when modeling the dorsal-side-only muscle inhibition, the parameter $H$ for describing max degree of optogenetic muscle inhibition was modulated to $H = 0.5 * H_{o p t i m a l}$ to qualitatively agree with experimental observations (Figure 6). This factor accounts for unequal degrees of paralysis during ventral vs. dorsal illumination (Figure 6—figure supplement 1), causing the PRC of dorsal-side illumination to show a relatively moderate response compared to ventral-side illumination.

To simulate the muscle inhibition on our threshold-switch model, we multiplied $M_{a}$ with $(1 - Q)$ any time the model was to be inhibited during its oscillatory period. To apply this operation to the alternative models, we factored out a term as a generalized active muscle moment for each individual model and then multiplied it with the bell-shaped function described above. The generalized forms of active muscle moment for the alternative models are implemented by modifying their original forms as follows:

a. For the van der Pol Oscillator, it is modified as:

{\begin{cases} \ddot{K} + (- \tilde{M_{V}} + P_{V}) \dot{K} + ω_{V}^{2} K = 0 \\ M_{V} = a_{V} (1 - b_{V} K^{2}) + P_{V} \end{cases};

b. For the Rayleigh Oscillator, it is modified as:

{\begin{cases} \ddot{K} + (- \tilde{M_{R}} + P_{R}) \dot{K} + ω_{R}^{2} K = 0 \\ M_{R} = a_{R} (1 - b_{R} {\dot{K}}^{2}) + P_{R} \end{cases};

c. For the Stuart-Landau Oscillator, it is modified as:

{\begin{cases} \dot{Z} + (- \tilde{M_{S}} + P_{S}) Z = 0 \\ M_{S} = σ - \frac{l}{2} {| Z |}^{2} + P_{S} \end{cases} .

For each individual model listed above, $\tilde{M_{i}}$ (subscript $i$ represents V, R, and S, respectively) is the generalized muscle moment which is to be multiplied by the bell-shaped factor $(1 - Q)$ upon perturbation, and $P_{i}$ is the additional damping coefficient. Note that the minus sign prior to $M_{i}$ in the first equation of each set indicates that $M_{i}$ is a negative damping term that provides power to the system, while $P_{i}$ is set positive for modeling the effect of bending toward the straight posture due to internal and external viscosity. Also note that Equations A15-A17 would be equivalent to their original form (Equations A11–A13) when inhibition is absent (in this case, $\tilde{M_{i}} = M_{i}$ ).

By modeling the muscle inhibition process during locomotion, we were able to perform simulations of phase response experiments on individual models to produce perturbed systems dynamics (Figure 5—figure supplement 1J–L) and the corresponding PRCs (Figure 5—figure supplement 1N–P and Figure 6—figure supplement 2).

Optimization of models

For each individual model we developed, the parameters were determined via a two-round fitting process. First, a subset of parameters was determined by fitting the model to observations of free-moving dynamics; the model could generate free-moving dynamics close to observations at this point. Second, the rest of the parameters were settled by fitting it to experimental phase response curves; a model would be fully determined at this point. Detailed descriptions of the two-step optimization procedure for individual models are provided as follows:

For the original threshold-switch model, parameters $τ_{u}$ , $M_{0}$ , and $τ_{m}$ were explicitly estimated from the experiments of free locomotion using phase portrait techniques described above. To simulate free locomotion, we further determined the position of switch points in the model (as indicated in Figure 4C red circle), which we did using method described by Equation A7. Next, we plugged the determined parameters into the model and conducted the second round of optimization by fitting the model with undetermined parameters $P_{t h}$ , $b$ , as well as the parameters for simulating muscle inhibition — $H$ and $q$ . We generated model PRC by perturbing the model oscillator at different times within a cycle period and settled the parameters such that the model PRC matched the experimental one with a minimum mean squared error (MSE) (During the computation of MSE, values of both model and experimental PRCs were sampled across the entire range of $ϕ$ with 100 evenly distributed samples. In this case, $Δ ϕ = 2 π / 100$ ):

\begin{matrix} (P_{t h}, b, H, q) = \underset{P_{t h}, b, H, q}{a r g m i n} \sum_{0}^{2 π} {(P R C_{m o d e l} (P_{t h}, b, H, q; ϕ) - P R C_{e x p e r i m e n t} (ϕ))}^{2} Δ ϕ \end{matrix}

To find the parameters that minimize the difference, a global minimum search was performed using the MATLAB function ‘GlobalSearch’ (Ugray et al., 2007). When run, the function repeatedly uses a local minimum solver with different batches of parameter range and attempts to locate a solution that produces the lowest MSE value.

Similarly, the two-step optimization procedures for individual alternative models are summarized in Appendix 1—table 1.

Appendix 1—table 1

Objective functions used in the optimization procedures for alternative models.

Type	Free locomotion model	Complete model
van der Pol	$\underset{a_{V}, b_{V}, ω_{V}}{a r g m i n} ({(\frac{T_{v d P}}{T_{e x p t}} - 1)}^{2} + {(\frac{A_{v d P}}{A_{e x p}} - 1)}^{2})$	$\underset{p_{V}, H, q}{a r g m i n} \sum_{0}^{2 π} (P R C_{v d p} (ϕ) - P R C_{e x p} (ϕ))^{2} Δ ϕ$
Rayleigh	$\underset{a_{R}, b_{R}, ω_{R}}{a r g m i n} ({(\frac{T_{R a y l e i g h}}{T_{e x p t}} - 1)}^{2} + {(\frac{A_{R a y l e i g h}}{A_{e x p}} - 1)}^{2})$	$\underset{p_{R}, H, q}{a r g m i n} \sum_{0}^{2 π} (P R C_{R a y l e i g h} (ϕ) - P R C_{e x p} (ϕ))^{2} Δ ϕ$
Stuart-Landau	$\underset{a_{S}, b_{S}, ω_{S}}{a r g m i n} ({(\frac{T_{S L}}{T_{e x p t}} - 1)}^{2} + {(\frac{A_{S L}}{A_{e x p}} - 1)}^{2})$	$\underset{p_{S}, H, q}{a r g m i n} \sum_{0}^{2 π} (P R C_{S L} (ϕ) - P R C_{e x p} (ϕ))^{2} Δ ϕ$

Two-step optimization procedure for van der Pol, Rayleigh, and Stuart-Landau oscillators. The first-step optimization determines part of parameters such that individual models generate free locomotion dynamics. The second-step optimization leads to complete models such that models’ perturbed dynamics and phase response curves are produced.

Data availability

All data and software have been deposited to a Dryad repository (https://doi.org/10.5061/dryad.wwpzgmsk2).

The following data sets were generated

1. Ji H
2. Fouad AD
3. Teng S
4. Liu A
5. Alvarez-Illera P
6. Yao B
7. Li Z
8. Fang-Yen C
(2021) Dryad Digital Repository
Data from: Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in Caenorhabditis elegans.

https://doi.org/10.5061/dryad.wwpzgmsk2

References

(2005) The Kuramoto model: a simple paradigm for synchronization phenomena
Reviews of Modern Physics 77:137–185.

https://doi.org/10.1103/RevModPhys.77.137
- Google Scholar
1. Ache JM
2. Matheson T
(2013) Passive joint forces are tuned to limb use in insects and drive movements without motor activity
Current Biology 23:1418–1426.

https://doi.org/10.1016/j.cub.2013.06.024
- PubMed
- Google Scholar
1. Akitake B
2. Ren Q
3. Boiko N
4. Ni J
5. Sokabe T
6. Stockand JD
7. Eaton BA
8. Montell C
(2015) Coordination and fine motor control depend on Drosophila TRPγ
Nature Communications 6:1–13.

https://doi.org/10.1038/ncomms8288
- Google Scholar
(1981) Peripheral feedback mechanisms acting on the central pattern generators for locomotion in fish and cat
Canadian Journal of Physiology and Pharmacology 59:713–726.

https://doi.org/10.1139/y81-108
- PubMed
- Google Scholar
(1999) The Caenorhabditis elegans unc-49 locus encodes multiple subunits of a heteromultimeric GABA receptor
The Journal of Neuroscience 19:5348–5359.

https://doi.org/10.1523/JNEUROSCI.19-13-05348.1999
- PubMed
- Google Scholar
1. Bargmann CI
(1998) Neurobiology of the Caenorhabditis elegans genome
Science 282:2028–2033.

https://doi.org/10.1126/science.282.5396.2028
- PubMed
- Google Scholar
1. Bässler U
(1977) Sensory control of leg movement in the stick insect Carausius morosus
Biological Cybernetics 25:61–72.

https://doi.org/10.1007/BF00337264
- PubMed
- Google Scholar
1. Berri S
2. Boyle JH
3. Tassieri M
4. Hope IA
5. Cohen N
(2009) Forward locomotion of the nematode C. elegans is achieved through modulation of a single gait
HFSP Journal 3:186–193.

https://doi.org/10.2976/1.3082260
- PubMed
- Google Scholar
(2009) Sensory feedback induced by front-leg stepping entrains the activity of central pattern generators in caudal segments of the stick insect walking system
Journal of Neuroscience 29:2972–2983.

https://doi.org/10.1523/JNEUROSCI.3155-08.2009
- PubMed
- Google Scholar
1. Boyle JH
2. Berri S
3. Cohen N
(2012) Gait modulation in C. elegans: an integrated neuromechanical model
Frontiers in Computational Neuroscience 6:10.

https://doi.org/10.3389/fncom.2012.00010
- PubMed
- Google Scholar
1. Brenner S
(1974) The genetics of Caenorhabditis elegans
Genetics 77:71–94.

https://doi.org/10.1093/genetics/77.1.71
- PubMed
- Google Scholar
1. Brodfuehrer PD
2. Friesen WO
(1986) From stimulation to undulation: a neuronal pathway for the control of swimming in the leech
Science 234:1002–1004.

https://doi.org/10.1126/science.234.4779.1002
- PubMed
- Google Scholar
1. Bryden J
2. Cohen N
(2008) Neural control of Caenorhabditis elegans forward locomotion: the role of sensory feedback
Biological Cybernetics 98:339–351.

https://doi.org/10.1007/s00422-008-0212-6
- PubMed
- Google Scholar
1. Butler VJ
2. Branicky R
3. Yemini E
4. Liewald JF
5. Gottschalk A
6. Kerr RA
7. Chklovskii DB
8. Schafer WR
(2015) A consistent muscle activation strategy underlies crawling and swimming in Caenorhabditis elegans
Journal of the Royal Society Interface 12:20140963.

https://doi.org/10.1098/rsif.2014.0963
- Google Scholar
1. Cang J
2. Yu X
3. Friesen OW
(2001) Sensory modification of leech swimming: interactions between ventral stretch receptors and swim-related neurons
Journal of Comparative Physiology A 187:569–579.

https://doi.org/10.1007/s003590100229
- Google Scholar
1. Cang J
2. Friesen WO
(2000) Sensory modification of leech swimming: rhythmic activity of ventral stretch receptors can change intersegmental phase relationships
The Journal of Neuroscience 20:7822–7829.

https://doi.org/10.1523/JNEUROSCI.20-20-07822.2000
- PubMed
- Google Scholar
1. Cestnik R
2. Rosenblum M
(2018) Inferring the phase response curve from observation of a continuously perturbed oscillator
Scientific Reports 8:1–10.

https://doi.org/10.1038/s41598-018-32069-y
- PubMed
- Google Scholar
(1985) The neural circuit for touch sensitivity in Caenorhabditis elegans
The Journal of Neuroscience 5:956–964.

https://doi.org/10.1523/JNEUROSCI.05-04-00956.1985
- PubMed
- Google Scholar
Book
(1988)
The Nematode Caenorhabditis elegans

New York: Cold Spring Harbor Laboratory Press.
- Google Scholar
(2011) Mechanisms underlying rhythmic locomotion: body–fluid interaction in undulatory swimming
Journal of Experimental Biology 214:561–574.

https://doi.org/10.1242/jeb.048751
- Google Scholar
1. Cohen AH
2. Wallen P
(1980) The neuronal correlate of locomotion in fish
Experimental Brain Research 41:11–18.

https://doi.org/10.1007/BF00236674
- Google Scholar
1. Cook SJ
2. Jarrell TA
3. Brittin CA
4. Wang Y
5. Bloniarz AE
6. Yakovlev MA
7. Nguyen KCQ
8. Tang LT
9. Bayer EA
10. Duerr JS
11. Bülow HE
12. Hobert O
13. Hall DH
14. Emmons SW
(2019) Whole-animal connectomes of both Caenorhabditis elegans sexes
Nature 571:63–71.

https://doi.org/10.1038/s41586-019-1352-7
- PubMed
- Google Scholar
1. Croll NA
(1971) The behaviour of Nematodes, their activity, senses and responses
The Journal of Applied Ecolog 8:956.

https://doi.org/10.2307/2402696
- Google Scholar
1. Cruse H
(1976) The control of body position in the stick insect (Carausius morosus), when walking over uneven surfaces
Biological Cybernetics 24:25–33.

https://doi.org/10.1007/BF00365591
- Google Scholar
1. Deng L
2. Denham JE
3. Arya C
4. Yuval O
5. Cohen N
6. Haspel G
(2021) Inhibition underlies fast undulatory locomotion in Caenorhabditis elegans
Eneuro 8:ENEURO.0241-20.2020.

https://doi.org/10.1523/ENEURO.0241-20.2020
- PubMed
- Google Scholar
(2018) Signatures of proprioceptive control in Caenorhabditis elegans locomotion
Philosophical Transactions of the Royal Society B: Biological Sciences 373:20180208.

https://doi.org/10.1098/rstb.2018.0208
- Google Scholar
1. Der pol b VAN
(1940) Biological rhythms considered as relaxation oscillations
Acta Medica Scandinavica 103:76–88.

https://doi.org/10.1111/j.0954-6820.1940.tb11082.x
- Google Scholar
Book
1. Driscoll M
2. Kaplan J
(1997)
Mechanotransduction

In: Kaplan J, editors. C. elegans II. New York, United States: Cold Spring Harbor Laboratory Press. pp. 10–22.
- Google Scholar
1. Eastwood AL
2. Sanzeni A
3. Petzold BC
4. Park SJ
5. Vergassola M
6. Pruitt BL
7. Goodman MB
(2015) Tissue mechanics govern the rapidly adapting and symmetrical response to touch
PNAS 112:E6955–E6963.

https://doi.org/10.1073/pnas.1514138112
- PubMed
- Google Scholar
1. Ekeberg O
2. Grillner S
(1999) Simulations of neuromuscular control in lamprey swimming
Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 354:895–902.

https://doi.org/10.1098/rstb.1999.0441
- PubMed
- Google Scholar
1. Fang-Yen C
2. Wyart M
3. Xie J
4. Kawai R
5. Kodger T
6. Chen S
7. Wen Q
8. Samuel AD
(2010) Biomechanical analysis of gait adaptation in the nematode Caenorhabditis elegans
PNAS 107:20323–20328.

https://doi.org/10.1073/pnas.1003016107
- PubMed
- Google Scholar
Book
(1993) Statistical Analysis of Spherical Data
Cambridge University Press.

https://doi.org/10.1017/CBO9780511623059
- Google Scholar
1. Fitzhugh R
(1961) Impulses and physiological states in theoretical models of nerve membrane
Biophysical Journal 1:445–466.

https://doi.org/10.1016/S0006-3495(61)86902-6
- PubMed
- Google Scholar
1. Fouad AD
2. Teng S
3. Mark JR
4. Liu A
5. Alvarez-Illera P
6. Ji H
7. Du A
8. Bhirgoo PD
9. Cornblath E
10. Guan SA
11. Fang-Yen C
(2018) Distributed rhythm generators underlie Caenorhabditis elegans forward locomotion
eLife 7:e29913.

https://doi.org/10.7554/eLife.29913
- PubMed
- Google Scholar
Book
1. Friesen WO
(2009)
Central pattern generators: sensory feedback

In: Friesen WO, editors. Encyclopedia of Neuroscience. Elsevier. pp. 701–709.
- Google Scholar
1. Gao S
2. Guan SA
3. Fouad AD
4. Meng J
5. Kawano T
6. Huang YC
7. Li Y
8. Alcaire S
9. Hung W
10. Lu Y
11. Qi YB
12. Jin Y
13. Alkema M
14. Fang-Yen C
15. Zhen M
(2018) Excitatory motor neurons are local oscillators for backward locomotion
eLife 7:e29915.

https://doi.org/10.7554/eLife.29915
- PubMed
- Google Scholar
1. Graham D
(1985)
Advances in Insect Physiology

31–140, Pattern and control of walking in insects, Advances in Insect Physiology, Elsevier, 10.1016/S0065-2806(08)60039-9.
- Google Scholar
1. Gray J
(1933) Studies in animal locomotion: I. the movement of fish with special reference to the eel
Journal of Experimental Biology 10:88–104.

https://doi.org/10.1242/jeb.10.1.88
- Google Scholar
1. Gray J
2. Lissmann HW
(1964) The locomotion of nematodes
Journal of Experimental Biology 41:135–154.

https://doi.org/10.1242/jeb.41.1.135
- Google Scholar
(1981) Entrainment of the spinal pattern generators for swimming by mechano-sensitive elements in the lamprey spinal cord in vitro
Brain Research 217:380–386.

https://doi.org/10.1016/0006-8993(81)90015-9
- PubMed
- Google Scholar
1. Grillner S
(2003) The motor infrastructure: from ion channels to neuronal networks
Nature Reviews Neuroscience 4:573–586.

https://doi.org/10.1038/nrn1137
- PubMed
- Google Scholar
1. Grillner S
2. Wallén P
(2002) Cellular bases of a vertebrate locomotor system-steering, intersegmental and segmental co-ordination and sensory control
Brain Research Reviews 40:92–106.

https://doi.org/10.1016/S0165-0173(02)00193-5
- PubMed
- Google Scholar
1. Guo ZV
2. Mahadevan L
(2008) Limbless undulatory propulsion on land
PNAS 105:3179–3184.

https://doi.org/10.1073/pnas.0705442105
- PubMed
- Google Scholar
(2005) Phase-response curves give the responses of neurons to transient inputs
Journal of Neurophysiology 94:1623–1635.

https://doi.org/10.1152/jn.00359.2004
- PubMed
- Google Scholar
1. Hobert O
(2003) Behavioral plasticity in C. elegans: paradigms, circuits, genes
Journal of Neurobiology 54:203–223.

https://doi.org/10.1002/neu.10168
- PubMed
- Google Scholar
Book
1. Izhikevich EM
(2007)
Dynamical Systems in Neuroscience

Massachusetts, United States: MIT Press.
- Google Scholar
1. Izquierdo EJ
2. Beer RD
(2018) From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegans
Philosophical Transactions of the Royal Society B: Biological Sciences 373:20170374.

https://doi.org/10.1098/rstb.2017.0374
- Google Scholar
(2021) Neuromechanical Mechanisms of Gait Adaptation in C. elegans : Relative Roles of Neural and Mechanical Coupling
SIAM Journal on Applied Dynamical Systems 20:1022–1052.

https://doi.org/10.1137/20M1346122
- Google Scholar
(2020) Nested neuronal dynamics orchestrate a behavioral hierarchy across timescales
Neuron 105:562–576.

https://doi.org/10.1016/j.neuron.2019.10.037
- PubMed
- Google Scholar
(2008) Systems level circuit model of C. elegans undulatory locomotion: mathematical modeling and molecular genetics
Journal of Computational Neuroscience 24:253–276.

https://doi.org/10.1007/s10827-007-0054-6
- PubMed
- Google Scholar
(1975) On the movement of slender bodies near plane boundaries at low Reynolds number
Journal of Fluid Mechanics 72:529–540.

https://doi.org/10.1017/S0022112075003126
- Google Scholar
1. Kawano T
2. Po MD
3. Gao S
4. Leung G
5. Ryu WS
6. Zhen M
(2011) An imbalancing act: gap junctions reduce the backward motor circuit activity to Bias C. elegans for forward locomotion
Neuron 72:572–586.

https://doi.org/10.1016/j.neuron.2011.09.005
- PubMed
- Google Scholar
Book
1. Kiehn O
(1998)
Neuronal Mechanisms for Generating Locomotor Activity

New York, United States: New York Academy of Sciences.
- Google Scholar
1. Kiehn O
(2011) Development and functional organization of spinal locomotor circuits
Current Opinion in Neurobiology 21:100–109.

https://doi.org/10.1016/j.conb.2010.09.004
- Google Scholar
1. Kristan WB
2. Calabrese RL
(1976) Rhythmic swimming activity in neurones of the isolated nerve cord of the leech
Journal of Experimental Biology 65:643–668.

https://doi.org/10.1242/jeb.65.3.643
- Google Scholar
(2017) Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotion
PLOS Computational Biology 13:e1005303.

https://doi.org/10.1371/journal.pcbi.1005303
- PubMed
- Google Scholar
(2011) Optogenetic manipulation of neural activity in freely moving Caenorhabditis elegans
Nature Methods 8:147–152.

https://doi.org/10.1038/nmeth.1554
- PubMed
- Google Scholar
1. Li W
2. Feng Z
3. Sternberg PW
4. Xu XZ
(2006) A C. elegans stretch receptor neuron revealed by a mechanosensitive TRP channel homologue
Nature 440:684–687.

https://doi.org/10.1038/nature04538
- PubMed
- Google Scholar
1. Marder E
2. Bucher D
(2001) Central pattern generators and the control of rhythmic movements
Current Biology 11:R986–R996.

https://doi.org/10.1016/S0960-9822(01)00581-4
- PubMed
- Google Scholar
1. Marder E
2. Calabrese RL
(1996) Principles of rhythmic motor pattern generation
Physiological Reviews 76:687–717.

https://doi.org/10.1152/physrev.1996.76.3.687
- PubMed
- Google Scholar
(1993) The GABAergic nervous system of Caenorhabditis elegans
Nature 364:337–341.

https://doi.org/10.1038/364337a0
- PubMed
- Google Scholar
1. Mendes CS
2. Bartos I
3. Akay T
4. Márka S
5. Mann RS
(2013) Quantification of gait parameters in freely walking wild type and sensory deprived Drosophila melanogaster
eLife 2:e00231.

https://doi.org/10.7554/eLife.00231
- PubMed
- Google Scholar
(1997) Contractile properties of obliquely striated muscle from the mantle of squid (Alloteuthis subulata) and cuttlefish (Sepia officinalis)
Journal of Experimental Biology 200:2425–2436.

https://doi.org/10.1242/jeb.200.18.2425
- PubMed
- Google Scholar
1. Morris C
2. Lecar H
(1981) Voltage oscillations in the barnacle giant muscle fiber
Biophysical Journal 35:193–213.

https://doi.org/10.1016/S0006-3495(81)84782-0
- PubMed
- Google Scholar
(2011) Neuronal control of swimming behavior: comparison of vertebrate and invertebrate model systems
Progress in Neurobiology 93:244–269.

https://doi.org/10.1016/j.pneurobio.2010.11.001
- PubMed
- Google Scholar
Book
1. Nadim F
2. Zhao S
3. Bose A
(2012) A PRC Description of How Inhibitory Feedback Promotes Oscillation stability
In: Bose A, editors. Phase Response Curves in Neuroscience. Springer. pp. 399–417.

https://doi.org/10.1007/978-1-4614-0739-3_16
- Google Scholar
(1962) An active pulse transmission line simulating nerve axon
Proceedings of the IRE 50:2061–2070.

https://doi.org/10.1109/JRPROC.1962.288235
- Google Scholar
Report
1. Nave CR
(2007) Hyperphysics
Georgia State University.

http://hyperphysics.phy-astr.gsu.edu/hbase/hph.html
- Google Scholar
1. Niebur E
2. Erdös P
(1991) Theory of the locomotion of Nematodes: dynamics of undulatory progression on a surface
Biophysical Journal 60:1132–1146.

https://doi.org/10.1016/S0006-3495(91)82149-X
- PubMed
- Google Scholar
(2005) The MEC-4 DEG/ENaC channel of Caenorhabditis elegans touch receptor neurons transduces mechanical signals
Nature Neuroscience 8:43–50.

https://doi.org/10.1038/nn1362
- PubMed
- Google Scholar
(2021) A neuromechanical model of multiple network rhythmic pattern generators for forward locomotion in C. elegans
Frontiers in Computational Neuroscience 15:10.

https://doi.org/10.3389/fncom.2021.572339
- Google Scholar
1. Pearce RA
2. Friesen WO
(1984) Intersegmental coordination of leech swimming: comparison of in situ and isolated nerve cord activity with body wall movement
Brain Research 299:363–366.

https://doi.org/10.1016/0006-8993(84)90720-0
- PubMed
- Google Scholar
Book
1. Pearson KG
(2004) Generating the Walking Gait: Role of Sensory feedbackProgress in Brain Research
Amsterdam, Netherlands: Elsevier.

https://doi.org/10.1016/S0079-6123(03)43012-4
- Google Scholar
1. Pierce-Shimomura JT
2. Chen BL
3. Mun JJ
4. Ho R
5. Sarkis R
6. McIntire SL
(2008) Genetic analysis of crawling and swimming locomotory patterns in C. elegans
PNAS 105:20982–20987.

https://doi.org/10.1073/pnas.0810359105
- PubMed
- Google Scholar
1. Pietras B
2. Daffertshofer A
(2019) Network dynamics of coupled oscillators and phase reduction techniques
Physics Reports 819:1–105.

https://doi.org/10.1016/j.physrep.2019.06.001
- Google Scholar
Book
1. Rayleigh J
(1896)
The Theory of Sound

London: Macmillan.
- Google Scholar
1. Roberts TJ
2. Azizi E
(2011) Flexible mechanisms: the diverse roles of biological springs in vertebrate movement
Journal of Experimental Biology 214:353–361.

https://doi.org/10.1242/jeb.038588
- Google Scholar
(2010) Kinetic properties of mechanically activated currents in spinal sensory neurons
The Journal of Physiology 588:301–314.

https://doi.org/10.1113/jphysiol.2009.182360
- PubMed
- Google Scholar
Book
(2011) Phase Response Curves in Neuroscience
Berlin, Germany: Springer.

https://doi.org/10.1007/978-1-4614-0739-3
- Google Scholar
(2010) Phase-response curves and synchronized neural networks
Philosophical Transactions of the Royal Society B: Biological Sciences 365:2407–2422.

https://doi.org/10.1098/rstb.2009.0292
- Google Scholar
1. Stiernagle T
(2006) Maintenance of C. elegans
WormBook 11:1–11.

https://doi.org/10.1895/wormbook.1.101.1
- Google Scholar
Book
1. Sulston J
2. Hodgkin J
(1988)
The Nematode Caenorhabditis Elegans

CSHL Press.
- Google Scholar
1. Tao L
2. Porto D
3. Li Z
4. Fechner S
5. Lee SA
6. Goodman MB
7. Xu XZS
8. Lu H
9. Shen K
(2019) Parallel processing of two mechanosensory modalities by a single neuron in C. elegans
Developmental Cell 51:617–631.

https://doi.org/10.1016/j.devcel.2019.10.008
- PubMed
- Google Scholar
1. Ugray Z
2. Lasdon L
3. Plummer J
4. Glover F
5. Kelly J
6. Martí R
(2007) Scatter search and local NLP solvers: a multistart framework for global optimization
INFORMS Journal on Computing 19:328–340.

https://doi.org/10.1287/ijoc.1060.0175
- Google Scholar
1. Van der Pol B
(1926) Lxxxviii.on “relaxation-oscillations”
The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2:978–992.

https://doi.org/10.1080/14786442608564127
- Google Scholar
(2006) The motor circuit
International Review of Neurobiology 69:125–167.

https://doi.org/10.1016/S0074-7742(05)69005-8
- PubMed
- Google Scholar
1. Wallace HR
(1968) The dynamics of nematode movement
Annual Review of Phytopathology 6:91–114.

https://doi.org/10.1146/annurev.py.06.090168.000515
- Google Scholar
1. Wen Q
2. Po MD
3. Hulme E
4. Chen S
5. Liu X
6. Kwok SW
7. Gershow M
8. Leifer AM
9. Butler V
10. Fang-Yen C
11. Kawano T
12. Schafer WR
13. Whitesides G
14. Wyart M
15. Chklovskii DB
16. Zhen M
17. Samuel AD
(2012) Proprioceptive coupling within motor neurons drives C. elegans forward locomotion
Neuron 76:750–761.

https://doi.org/10.1016/j.neuron.2012.08.039
- PubMed
- Google Scholar
1. Wendler G
(1968)
Ein Analogmodell Der Beinbewegungen Eines laufenden insekts

Kybernetik 18:68–74.
- Google Scholar
(1986) The structure of the nervous system of the nematode Caenorhabditis elegans
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 314:1–340.

https://doi.org/10.1098/rstb.1986.0056
- PubMed
- Google Scholar
Book
1. Winfree AT
(2001)
The Geometry of Biological Time

Springer Science & Business Media.
- Google Scholar
1. Xu T
2. Huo J
3. Shao S
4. Po M
5. Kawano T
6. Lu Y
7. Wu M
8. Zhen M
9. Wen Q
(2018) Descending pathway facilitates undulatory wave propagation in Caenorhabditis elegans through gap junctions
PNAS 115:E4493–E4502.

https://doi.org/10.1073/pnas.1717022115
- PubMed
- Google Scholar
1. Yeon J
2. Kim J
3. Kim DY
4. Kim H
5. Kim J
6. Du EJ
7. Kang K
8. Lim HH
9. Moon D
10. Kim K
(2018) A sensory-motor neuron type mediates proprioceptive coordination of steering in C. elegans via two TRPC channels
PLOS Biology 16:e2004929.

https://doi.org/10.1371/journal.pbio.2004929
- PubMed
- Google Scholar
1. Yu X
2. Nguyen B
3. Friesen WO
(1999) Sensory feedback can coordinate the swimming activity of the leech
The Journal of Neuroscience 19:4634–4643.

https://doi.org/10.1523/JNEUROSCI.19-11-04634.1999
- PubMed
- Google Scholar
1. Yu X
2. Friesen WO
(2004) Entrainment of leech swimming activity by the ventral stretch receptor
Journal of Comparative Physiology A 190:939–949.

https://doi.org/10.1007/s00359-004-0549-9
- Google Scholar

Article and author information

Author details

Hongfei Ji

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

Contributed equally with
Anthony D Fouad

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9617-6411
Anthony D Fouad

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Software, Investigation, Methodology, Writing - original draft, Writing - review and editing

Contributed equally with
Hongfei Ji

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4677-2968
Shelly Teng

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Investigation

Competing interests
No competing interests declared
Alice Liu

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Investigation

Competing interests
No competing interests declared
Pilar Alvarez-Illera

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Investigation

Competing interests
No competing interests declared
Bowen Yao

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Software, Investigation

Competing interests
No competing interests declared
Zihao Li

Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States

Contribution
Investigation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4304-8322
Christopher Fang-Yen
1. Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania, Philadelphia, United States
2. Department of Neuroscience, Perelman School of Medicine, University of Pennsylvania, Philadelphia, United States
Contribution
Conceptualization, Resources, Software, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing

For correspondence
fangyen@seas.upenn.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4568-3218

Funding

National Institutes of Health (1R01NS084835)

Hongfei Ji
Anthony D Fouad
Christopher Fang-Yen

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Mei Zhen and Quan Wen for providing strains. Some strains were provided by the CGC, funded by NIH Office of Research Infrastructure Programs (P40 OD010440). We thank Gal Haspel, Michael Carchidi, and Patrick McClanahan for helpful discussions. HJ, ADF, and CF-Y. were supported by the National Institutes of Health (1R01NS084835). ST was supported by an Abraham Noordergraaf Research Fellowship and a Littlejohn Fellowship.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.