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

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 cord . 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.

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 k 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 ¼ k Á 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).
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 p=4 radians relative to the curvature, and the resulting phase portrait would be symmetric about both the K and dK=dt 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 dK=dt axes ( Figure 2D). Plots of both the phase portrait ( Figure 2D) and the time dependence ( Figure 2C) show that K and dK=dt 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 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).
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; p=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 p=3; p=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 2p=3; 5p=6 ½ , the head response was similar to that within the interval 0; p=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 ( 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 f » p=3 and 4p=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, Figure 3 continued (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 p=3; p=2 ½ . (I-K) Similar to (C-E), for phase range 2p=3; 5p=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:16p 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. The online version of this article includes the following figure supplement(s) for figure 3:        Video 1. Transient illumination of the anterior region of a freely moving Pmyo-3::NpHR worm. Green-shaded region indicates timing and location of illumination.
https://elifesciences.org/articles/69905#video1 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 t m » 100 ms. This time constant is much smaller than the duration of each muscle moment plateau period ( » 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: where t 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 (dK=dt ¼ 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 th and ÀP th ( Figure 4D), the active muscle moment undergoes a change of sign ( Figure 4E), causing the head to bend toward the opposite direction ( Figure 4B).
Our model has 5 parameters: (1) t u , the bending relaxation time scale, (2) t m , the muscle switching time scale, (3) M 0 , the amplitude of the scaled active muscle moment, (4-5) b and P th , 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 th 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 t u , M 0 , and t 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 th , 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 th , the optimization estimated their values to be b ¼ 0:046 s and P th ¼ 2:33, as shown on the phase portraits (gray dashed lines in Figures 4C, 5B and D).
The threshold-switch mechanism model provides an explanation for the observed sawtoothshaped 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 singleside 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. 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. We incorporated the effect of external viscosity into our model through the bending relaxation time constant t u (see Appendix). We ran our model to determine the dependence of model output on viscosity with varying viscosity h. 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 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 th and ÀP th . 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 _ 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 wellknown 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, velocitydependent 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 5figure 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 meansquare 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 t u and the muscle moment transition time scale t 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 » 282ms, which is quite close to our measured result, t u » 260ms. Furthermore, our measurement of the muscle moment transition time scale (t m » 100ms) 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 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 curvaturesensing 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.

Worm strains and cultivation
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 . All experiments were performed with young adult (< 1 day) hermaphrodites synchronized by hypochlorite bleaching.

Locomotion and phase response analyses
To perform quantitative recordings of worm behavior, we used a custom-built optogenetic targeting system as previously described Leifer et al., 2011). Analysis of images for worm's body posture was performed using a previously developed custom software . 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 2 ) 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
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). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.  ], WEN001 wenIs001 ].

Author contributions
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 mL of a 100 mM solution of ATR in ethanol to an overnight culture of 250 mL 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-mm-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  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 . The worm centerline of each image was smoothed using a cubic spline fit. We calculated curvature k 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 k 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 2p. The averaged curvature within one cycle was then computed by averaging all individual cycles in the Similarly, the averaged phase derivative of curvature within one cycle was calculated as dK=df

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; _ K À Á in the unit of curvature where ¼ 0:135s 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 _ 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 _ 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 Here, D f ð Þ is distance to the center of oscillation on the phase plane, which is set to the origin, such where f denotes the phase value of the current state estimated by the four-quadrant inverse tangent of the variable pair K; _ K À Á . In this expression, D C f ð Þ denotes the distance to the center of oscillation that is evaluated exactly on the normal cycle at phase f.
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 f C t ð Þ ¼ ! 0 Á t modT0 , where ! 0 ¼ 2p=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 (f C ¼ 0) to be when K reaches local maximum (or Þ) and hence f C ¼ p to be when K reaches local minimum (or x ¼ K min ; 0 ð Þ). In this way, we parameterized the normal cycle by defining a bijective map between phases and state points The map F x ð Þ ¼ f 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 f 0 ¼ F 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 _ 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 _ 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 Þ and x f ¼ p ð Þ (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 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 Àp. Here, f x C N À Á was computed twice using phase states x f ¼ 0 ð Þ [subscripted with u] and x f ¼ p ð Þ [subscripted with l] as references, respectively; (2) we calculated the induced phase shift PRC t illum ð Þ and the phase of the notch point is f 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 Þ that are defined above and took derivative of the trajectories with respect to time. Thus, by collecting all the phase states K; _ K À Á and their corresponding time derivatives dK=dt; d _ K À Á =dt À Á 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 f, 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 illum ; 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 TZ À2 , TZ À1 , TZ þ1 , and TZ þ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 TH À2 , TH À1 , TH þ1 , and TH þ2 , in the ascending order of time. (iv) With these measurements, cycle period T 0 was computed as T 0 ¼ TZ þ2 À TZ þ1 þ TZ À1 À TZ À2 þ TH þ2 À TH þ1 þ TH À1 À TH À2 ð Þ=4, so the angular frequency of undulation ! 0 ¼ 2p=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 f of each individual trial was computed as , and the corresponding phase shift F was computed as 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; 2p ½ for x and range Àp; p ½ 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:16p along the direction of f from 0 to 2p. 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 f 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 (

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: 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, k » y ss , we arrive at: Under the assumptions of small-amplitude undulations and a fixed wavelength l down the worm body, k can be considered as a travelling sinusoidal wave with a small deviation, k s; t ð Þ ¼ k 0 sin 2ps=l À !t ð Þ þ d, which leads to an approximation, k ssss » 2p=l ð Þ 4 k. Plugging this approximation into Equation A2 while keeping s fixed, after some rearrangement, one gets: In terms of the dimensionless curvature K ¼ k Á L and dimensionless muscle moment we can rewrite Equation A3 as: and we note that Equations A5 and A6 yield Equation 1. In Equation A6, both the wavelength l and the normal viscous drag coefficient C N vary with the fluid viscosity h (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: i. Proprioceptive feedback is sensed as a linear combination of the current curvature value and the current rate of change of curvature, P ¼ K þ b _ K (black curve in Figure 4D). ii. During movement of bending, this proprioceptive feedback is constantly compared with two threshold values P th and ÀP th (gray dashed bars in Figure 4D). iii. 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). iv. 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: 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 j t¼0 ¼ M 0 and Kj 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 , t u , and t m were estimated from data of free locomotion using phase portrait techniques described in the following subsection. Parameters b and P th were yet degenerate in this model of free locomotion. Here, we temporarily set b ¼ 0 and then set P th such that the oscillatory period predicted by model matched the average period measured from experiments with a minimum squared error: The nondegeneracy of b and P th 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, t u ¼ 260 ms, t m ¼ 100 ms, b ¼ 46 ms, and P th ¼ 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 þ t u _ 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 _ K to identify the two 'flat' regions and, through linear fits, obtained two linear relations respectively: . Thus, the bending relaxation time scale t u and the amplitude of the scaled active muscle moment are estimated ast 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 t u (estimated from the above method), K h i and _ 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): To estimate t 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: We notice that when t ¼ 0, the maximum of jdM a =dtj is achieved and the value is M 0 =t m . On the other hand, the maximum of dM a =dt j j can be obtained from the experimental observation by simply finding the peak of |dM a =dt| curve where M a ¼ K t ð Þ h i þ t u Á dK t ð Þ=dt h i. Thus, t m can be estimated as:

Parameter estimation
For our original threshold-switch model, parameters t u , t 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 th describe the proprioceptive feedback and the threshold-switch features in our model. Specifically, they characterize two threshold lines, K þ b _ K ¼ AEP th (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 th such that the oscillatory period defined by model matched the average period measured from the experiments.
To obtain the nondegeneracy of P th 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 th 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 t m ¼ 100 ms, M 0 ¼ 8:45, b ¼ 46 ms, P th ¼ 2:33. For the exogenous parameters, only the time constant of undulation, t u , varies according to external conditions. According to Equation A6, t 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 Á 10 À16 Nm 3 s; the bending modulus of worm body is measured as a ¼ 9:5 Á 10 À14 Nm 3 ; C N ¼ 31h is the coefficient of viscous drag for movement normal to the body (Katz et al., 1975), where h 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=L ¼ À0:158 log 10 h=h 0 ð Þ þ 1:5 for a continuous model realization in undulatory frequency and amplitude. Here, l is the wavelength and h 0 ¼ 1 mPa Á s. 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 Dt ð Þ, as a function of the time interval Dt 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 bellshaped 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 Dt ð Þ, which has a bell-shaped profile (Figure 4-figure supplement 1): 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 Á 10 À1=q so that Qj Dt¼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 optimal 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: b. For the Rayleigh Oscillator, it is modified as: c. For the Stuart-Landau Oscillator, it is modified as: For each individual model listed above,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,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

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 freemoving 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 t u , M 0 , and t 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 th , 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 f with 100 evenly distributed samples. In this case, Df ¼ 2p=100): 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- Two-step optimization procedure for van der Pol, Rayleigh, and Stuart-Landau oscillators. The firststep 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.