Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in Caenorhabditis elegans
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 halfcenter 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).
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), wellmapped 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 anteriortoposterior 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 Btype motor neurons are required for forward locomotion (Chalfie et al., 1985). The GABAergic Dtype motor neurons provide dorsoventral crossinhibition 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 Dtype 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 Btype 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 Btype 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, TRP4, 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 nlp12, 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, sawtoothshaped 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 relaxationoscillation process, with switching based on proprioceptive feedback, underlies the worm’s rhythmic dorsalventral 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 (FangYen 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 timedependent 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 timevarying curvature along the centerline of the body (FangYen et al., 2010; Leifer et al., 2011; PierceShimomura 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; FangYen 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 $\kappa $ 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=\kappa \xb7L$, 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 (FangYen 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 $\pi /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 nonellipsoidal 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 nonsinusoidal.
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 closedloop 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 Pmyo3::NpHR). Simultaneous muscle inhibition on both sides causes C. elegans to straighten due to internal elastic forces (FangYen 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 $\left[0,\pi /6\right]$, 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 $[\pi /3,\pi /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\pi /3,5\pi /6]$, the head response was similar to that within the interval $[0,\pi /6]$, and also resulted in a phase delay (Figure 3I–K).
Combining the data from all phases of inhibition yielded a sawtoothshaped 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 sawtoothshaped 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 2D representation we found that the phase shifts display a piecewise, linear increasing dependence on the phase of inhibition with two abrupt jumps occurring at $\varphi \approx \pi /3$ and $4\pi /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 (Punc17::NpHR::ECFP) or Btype motor neurons (Pacr5::ArchmCherry). In both strains, we again observed sawtoothshaped PRCs (Figure 3—figure supplements 6 and 7), with variations only in the magnitudes of phase shifts. These experiments show that the sawtoothshaped 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 Dtype motor neurons provide a dorsoventral reciprocal inhibition of opposing muscles during locomotion. We asked whether the Dtype 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 unc49(e407), a lossoffunction mutant of GABA_{A} receptor that is required by the Dtype 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 Dtype motor neurons are not necessary for the motor rhythm generator to show the sawtoothshaped PRC.
Sawtoothshaped 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 switchlike 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 (FangYen 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 (FangYen 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 ${\tau}_{m}\approx 100\phantom{\rule{thinmathspace}{0ex}}ms$. This time constant is much smaller than the duration of each muscle moment plateau period ($\approx 0.5\phantom{\rule{thinmathspace}{0ex}}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 ${\tau}_{u}$ describes the time scale of bending relaxation and ${M}_{a}\left(t\right)$ is the timevarying 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) ${\tau}_{u}$, the bending relaxation time scale, (2) ${\tau}_{m}$, the muscle switching time scale, (3) ${M}_{0}$, the amplitude of the scaled active muscle moment, (45) $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 tworound 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 bellshaped 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 ${\tau}_{u}$, ${M}_{0}$, and ${\tau}_{m}$ that had been explicitly estimated from freemoving experiments, we performed a tworound 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\phantom{\rule{thinmathspace}{0ex}}s$ and ${P}_{th}=2.33$, as shown on the phase portraits (gray dashed lines in Figures 4C, 5B and D).
The thresholdswitch 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 singleside 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 (Pmyo3::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 dorsalside or ventralside 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 dorsalside illumination shows a smaller paralytic response than that of ventralside 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 dorsalside 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; FangYen 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 loaddependent gait adaptation.
We incorporated the effect of external viscosity into our model through the bending relaxation time constant ${\tau}_{u}$ (see Appendix). We ran our model to determine the dependence of model output on viscosity with varying viscosity $\eta $. We found that model results for frequency and amplitude dependence on viscosity of the external medium are in quantitative agreement with previous experimental results (FangYen 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}_{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\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 wellknown mathematical descriptions of oscillators (van der Pol, Rayleigh, and StuartLandau oscillators) and compared them with our original thresholdswitch 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 secondorder differential equation for a harmonic oscillator with a nonlinear, displacementdependent damping term (see Appendix). By choosing a set of appropriate parameters, we found that the freerunning 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 bellshaped modulatory function over all phases within a cycle produced a similar sawtoothshaped 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 singleside muscle inhibitions to the system produced singlesawtoothshaped 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 selfsustained acoustic vibrations such as vibrating clarinet reeds (Rayleigh, 1896). It is based on a secondorder 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 freerunning 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 sawtoothshaped 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 StuartLandau 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 freerunning 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 freerunning 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 freerunning 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 thresholdswitch model (Figure 5—figure supplement 1M) or the van der Pol model (Figure 5—figure supplement 1N). Of all the models tested, the thresholdswitch model showed the least meansquare error with the PRC data (Figure 5—figure supplement 1M–P). We conclude that of these models, our thresholdswitch model produced the best overall agreement with experiments.
We also found that two important experimental findings, the nonsinusoidal freemoving dynamics and the sawtoothshaped PRCs can be achieved in our original model, the van der Pol and Rayleigh oscillators, which are all relaxation oscillators, but not in the StuartLandau 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 forwardmoving 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 crossinhibition 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 Btype motor neurons suffices to generate robust locomotory rhythms without crossinhibition.
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 ${\tau}_{u}$ and the muscle moment transition time scale ${\tau}_{m}$ (see Appendix for details), which may be compared with previous studies of worm biomechanics (FangYen et al., 2010; Berri et al., 2009) and neurophysiology (Milligan et al., 1997). FangYen 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 282ms$, which is quite close to our measured result, ${\tau}_{u}\approx 260ms$. Furthermore, our measurement of the muscle moment transition time scale (${\tau}_{m}\approx 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 subprocesses 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 selfsustained oscillatory system with intrinsic periodic relaxation/decay features (Van der Pol, 1926). The FitzhughNagumo 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 ‘pushpull’ 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 sawtoothshaped 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 Btype 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 Btype 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 rapidlyadapting (RA) and intermediateadapting (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 curvatureratesensing 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
Worm strains and cultivation
Request a detailed protocolC. 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 OP50ATR 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 protocolTo perform quantitative recordings of worm behavior, we used a custombuilt 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 highfrequency movements of the worm’s anterior tip. Descriptions of the apparatus and image analyses are available in Appendix.
For phase response experiments, opsinexpressing 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
Request a detailed protocolOur primary model is based on a novel neural control mechanism incorporated with a previously described biomechanical framework (FangYen 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[Pmyo3::NpHR::eCFP;lin15(+)];qhIs4[Pacr2::wCherry], YX119 unc49(e407);qhIs1[Pmyo3::NpHR::eCFP;lin15(+)], YX205 hpIs178[Punc17::NpHR::eCFP;lin15(+)], WEN001 wenIs001[Pacr5::Arch::mCherry;lin15(+)].
For optogenetic experiments, worms were cultivated in darkness on plates with OP50 containing the cofactor alltrans retinal (ATR). For control experiments and freemoving experiments, worms were cultivated on regular OP50 NGM plates without ATR. To make OP50ATR 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µmthick polyester shims (McMasterCarr 9513K42). For viscositydependence experiments, we used 5%, 17%, and 35% (by mass) solutions of dextran (SigmaAldrich 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 (FangYen 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 custombuilt 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 customwritten C++ software (Fouad et al., 2018) to perform realtime segmentation of the worm during image acquisition. The worm was identified in each image by its boundary and centerline calculated from a binary image. Anteriorposterior 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.
Postacquisition 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 $\kappa $ 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 $\kappa $ 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\left(t\right)$ was defined as the average of the dimensionless curvature over body coordinate 0.10.3; this range avoided highfrequency 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 timedependent to a phasedependent curvature by uniformly rescaling each cycle to a phase range of $2\pi $. The averaged curvature within one cycle was then computed by averaging all individual cycles in the phase domain: $\u27e8K\left(\varphi \right)\u27e9=\sum _{i=1}^{N}{K}_{i}\left(\varphi \right)/N$. Similarly, the averaged phase derivative of curvature within one cycle was calculated as $\u27e8dK/d\varphi \u27e9={\sum}_{i=1}^{N}(d{K}_{i}/d\varphi )/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 twodimensional variable, $\mathit{x}=\left(K,\xi \dot{K}\right)$ in the unit of curvature where $\xi =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 $\xi \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\left(\varphi \right)=\left(D\right(\varphi ){D}^{C}(\varphi \left)\right)/{D}^{C}\left(\varphi \right)$. Here, $D\left(\varphi \right)$ is distance to the center of oscillation on the phase plane, which is set to the origin, such that $D\left(\varphi \right)=\sqrt{K{\left(\varphi \right)}^{2}+{\left(\xi \dot{K}\left(\varphi \right)\right)}^{2}}$ where $\varphi $ denotes the phase value of the current state estimated by the fourquadrant inverse tangent of the variable pair $\left(K,\xi \dot{K}\right)$. In this expression, ${D}^{C}\left(\varphi \right)$ denotes the distance to the center of oscillation that is evaluated exactly on the normal cycle at phase $\varphi $.
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 ${\varphi}^{C}\left(t\right)={\omega}_{0}\xb7{t}_{mod{T}_{0}}$, where ${\omega}_{0}=2\pi /{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 (${\varphi}^{C}=0$) to be when $K$ reaches local maximum (or $\mathit{x}=({K}_{max},0)$) and hence ${\varphi}^{C}=\pi $ to be when $K$ reaches local minimum (or $\mathit{x}=({K}_{min},0)$). In this way, we parameterized the normal cycle by defining a bijective map between phases and state points $\mathrm{\Phi}\left({\mathit{x}}^{C}\right)={\varphi}^{C}$.
The map $\mathrm{\Phi}\left(\mathit{x}\right)=\varphi $ 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 ${\mathit{x}}_{0}$ is a point on the normal cycle and ${\mathit{y}}_{0}$ is a point off the normal cycle, then ${\mathit{y}}_{0}$ will have the same phase as ${\mathit{x}}_{0}$ if the trajectory starting at the initial point ${\mathit{y}}_{0}$ off the normal cycle converges to the trajectory starting at the initial point ${\mathit{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 ${\mathit{x}}_{0}$ on the normal cycle as the isochron (Winfree, 2001) for phase ${\varphi}_{0}=\mathrm{\Phi}\left({\mathit{x}}_{0}\right)$.
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' ${\mathit{x}}_{N}$ as it can be clearly seen from the phase plot (as shown in Figure 3E, H and K). After the notch point ${\mathit{x}}_{N}$, the oscillation will proceed to its next phase states $\mathit{x}\left(\varphi =0\right)$ and $\mathit{x}\left(\varphi =\pi \right)$ (or viceversa), both of which can be easily identified through peak finding from the curvature dynamics $K$. Hence, we obtained two subtrajectories from the oscillation: one ${\mathit{x}}_{N}\to \mathit{x}\left(\varphi =0\right)$, and the other ${\mathit{x}}_{N}\to \mathit{x}(\varphi =\pi )$. Next to determining the timing of the notch point $t\left({\mathit{x}}_{N}\right)$, 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\left({\mathit{x}}_{N}\right)$ as if the perturbation had not occurred, which is $\varphi}_{u}\left({\mathit{x}}_{N}^{C}\right)={\omega}_{0}{\left(t\left({\mathit{x}}_{N}\right)t\left(\mathit{x}\left(\varphi =0\right)\right)\right)}_{mod{T}_{0}$, or ${\varphi}_{l}\left({\mathit{x}}_{N}^{C}\right)={\omega}_{0}{\left(t\left({\mathit{x}}_{N}\right)t\left(\mathit{x}\left(\varphi =\pi \right)\right)\right)}_{mod{T}_{0}}\pi$. Here, $\varphi \left({\mathit{x}}_{N}^{C}\right)$ was computed twice using phase states $\mathit{x}\left(\varphi =0\right)$ [subscripted with $u$] and $\mathit{x}\left(\varphi =\pi \right)$ [subscripted with $l$] as references, respectively; (2) we calculated the induced phase shift $PRC\left({t}_{illum}\right)$ and the phase of the notch point is $\varphi \left({\mathit{x}}_{N}\right)=\varphi \left({\mathit{x}}_{N}^{C}\right)PRC\left({t}_{illum}\right)$. After obtaining the subtrajectories ${\mathit{x}}_{N}\to \mathit{x}\left(\varphi =0\right)$ and ${\mathit{x}}_{N}\to \mathit{x}(\varphi =\pi )$ and calculating the phase of ${\mathit{x}}_{N}$, we then estimated the phase values for all the points within each of the two subtrajectories 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 2D 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 2D 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 subtrajectories ${\mathit{x}}_{N}\to \mathit{x}\left(\varphi =0\right)$ and ${\mathit{x}}_{N}\to \mathit{x}(\varphi =\pi )$ that are defined above and took derivative of the trajectories with respect to time. Thus, by collecting all the phase states $\left(K,\xi \dot{K}\right)$ and their corresponding time derivatives $\left(dK/dt,d\left(\xi \dot{K}\right)/dt\right)$ 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 2D moving average for the raw outcome over the phase plane to smooth out the vector field. We used a linear 2D 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 $\varphi $, 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 timevarying 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 zerophase 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 halfcyclephase 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 ${\omega}_{0}=2\pi /{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 $\varphi $ of each individual trial was computed as ${\varphi}_{u}={\omega}_{0}{\left({T}_{illum}T{Z}_{1}\right)}_{mod{T}_{0}}$, ${\varphi}_{l}={\omega}_{0}{\left({T}_{illum}T{H}_{1}+{T}_{0}/2\right)}_{mod{T}_{0}}$, and the corresponding phase shift $F$ was computed as ${F}_{u}={\omega}_{0}{\left(T{Z}_{+1}T{Z}_{1}\right)}_{mod{T}_{0}}\pi $, ${F}_{l}={\omega}_{0}{\left(T{H}_{+1}T{H}_{1}+{T}_{0}/2\right)}_{mod{T}_{0}}\pi $. Here, phase of illumination and the corresponding phase shift were computed twice using zero [subscripted with $u$] and halfcycle [subscripted with $l$] phase points as references, respectively.
We generated 2D 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 2D bins with 25 bins on both dimensions covering the range $\left[0,2\pi \right]$ for x and range $\left[\pi ,\pi \right]$ for y. To indicate average tendency of phase shift depending on phase of illumination, we calculated a meancurve representation of PRCs via a moving average operation. In this process, each mean was calculated over a sliding window of width $0.16\pi $ along the direction of $\varphi $ from 0 to $2\pi $. 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 $\varphi $ 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 Pmyo3::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 sawtoothshape 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 (FangYen 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 (FangYen 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 smallamplitude approximation, $\kappa \approx {y}_{ss}$, we arrive at:
Under the assumptions of smallamplitude undulations and a fixed wavelength $\lambda $ down the worm body, $\kappa $ can be considered as a travelling sinusoidal wave with a small deviation, $\kappa \left(s,t\right)={\kappa}_{0}\mathrm{sin}\left(2\pi s/\lambda \omega t\right)+\delta $, which leads to an approximation, ${\kappa}_{ssss}\approx {\left(2\pi /\lambda \right)}^{4}\kappa $. Plugging this approximation into Equation A2 while keeping $s$ fixed, after some rearrangement, one gets:
In terms of the dimensionless curvature $K=\kappa \xb7L$ and dimensionless muscle moment
we can rewrite Equation A3 as:
where
and we note that Equations A5 and A6 yield Equation 1. In Equation A6, both the wavelength $\lambda $ and the normal viscous drag coefficient ${C}_{N}$ vary with the fluid viscosity $\eta $ (Berri et al., 2009 ; FangYen 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 relaxationoscillation 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}_{th}$ and ${P}_{th}$ (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 switchtriggered muscle transition we used a modified logistic function: ${M}_{a}\left(t\right)=\pm {M}_{0}\bullet \mathrm{tanh}\left(t/2{\tau}_{m}\right)$. Here, the plus sign indicates the dorsaltoventral 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\left(t\right)$ 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}$, ${\tau}_{u}$, and ${\tau}_{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$, ${\tau}_{u}=260\phantom{\rule{thinmathspace}{0ex}}ms$, ${\tau}_{m}=100\phantom{\rule{thinmathspace}{0ex}}ms$, $b=46\phantom{\rule{thinmathspace}{0ex}}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 freerunning 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+{\tau}_{u}\dot{K}={M}_{a}\left(t\right)$, the two flat regions indicate that the scaled active muscle moment, ${M}_{a}\left(t\right)$, 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: $\u27e8K\u27e9+{\tau}_{1}\u27e8\dot{K}\u27e9={M}_{1}$ (region 1) and $\u27e8K\u27e9+{\tau}_{2}\u27e8\dot{K}\u27e9={M}_{2}$ (region 2). Thus, the bending relaxation time scale ${\tau}_{u}$ and the amplitude of the scaled active muscle moment are estimated as $\hat{{\tau}_{u}}=\left({\tau}_{1}+{\tau}_{2}\right)/2$ and ${\hat{M}}_{\text{0}}=\left(\left{M}_{1}\right+\left{M}_{2}\right\right)/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 ${\tau}_{u}$ (estimated from the above method), $\u27e8K\u27e9$ and $\u27e8\dot{K}\u27e9$ (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}\left(t\right)$ on the right side of Equation 1. As expected and shown in Figure 2D Inset, the curve of ${M}_{a}\left(t\right)$ 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 ventraltodorsal (or viceversa) transition. We used a modified logistic function to model the ventraltodorsal muscle moment transition (substituting $t$ with $t$ for transition in the other direction):
To estimate ${\tau}_{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 $d{M}_{a}/dt$ is achieved and the value is ${M}_{0}/{\tau}_{m}$. On the other hand, the maximum of $\leftd{M}_{a}/dt\right$ can be obtained from the experimental observation by simply finding the peak of $d{M}_{a}/dt$ curve where ${M}_{a}=\u27e8K\left(t\right)\u27e9+\widehat{{\tau}_{u}}\cdot \u27e8dK\left(t\right)/dt\u27e9$. Thus, ${\tau}_{m}$ can be estimated as:
Parameter estimation
For our original thresholdswitch model, parameters ${\tau}_{u}$, ${\tau}_{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 thresholdswitch features in our model. Specifically, they characterize two threshold lines, $K+b\dot{K}=\pm {P}_{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 freemoving 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 ${\tau}_{m}=100\phantom{\rule{thinmathspace}{0ex}}ms$, ${M}_{0}=8.45$, $b=46\phantom{\rule{thinmathspace}{0ex}}ms$, ${P}_{th}=2.33$. For the exogenous parameters, only the time constant of undulation, ${\tau}_{u}$, varies according to external conditions. According to Equation A6, ${\tau}_{u}$ is explicitly determined in terms of other physical parameters, including biomechanical parameters measured in previous work (FangYen et al., 2010): the internal viscosity of worm body is measured as ${a}_{v}=5\cdot {10}^{16}\phantom{\rule{thinmathspace}{0ex}}N{m}^{3}s$; the bending modulus of worm body is measured as $a=9.5\cdot {10}^{14}\phantom{\rule{thinmathspace}{0ex}}N{m}^{3}$; ${C}_{N}=31\eta $ is the coefficient of viscous drag for movement normal to the body (Katz et al., 1975), where $\eta $ is the fluid viscosity. According to previous measurements of undulatory wavelengths in different viscous solutions (FangYen et al., 2010), we applied a logarithmic fit to the data points, yielding $\lambda /L=0.158{\mathrm{log}}_{10}(\eta /{\eta}_{0})+1.5$ for a continuous model realization in undulatory frequency and amplitude. Here, $\lambda $ is the wavelength and ${\eta}_{0}=1\phantom{\rule{thinmathspace}{0ex}}mPa\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 selfoscillator models described by 2D nonlinear systems: the van der Pol, Rayleigh, and StuartLandau oscillators.
First, we developed a model oscillator in the form taken from the van der Pol Oscillator:
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\phantom{\rule{thinmathspace}{0ex}}{s}^{1},{b}_{V}=2.04,{\omega}_{V}=5.51\phantom{\rule{thinmathspace}{0ex}}{s}^{1}.$
Second, we developed a model oscillator by modifying the Rayleigh Oscillator:
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\phantom{\rule{thinmathspace}{0ex}}{s}^{1},{b}_{R}=0.0023\phantom{\rule{thinmathspace}{0ex}}{s}^{2},{\omega}_{R}=5.6\phantom{\rule{thinmathspace}{0ex}}{s}^{1}.$
Third, we developed a model oscillator by modifying the StuartLandau Oscillator:
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 $\sigma $ 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 $\leftZ\right$. Simulated oscillation is a limit cycle of the model (Figure 5—figure supplement 1D,H), given parameters ${l}_{r}=0.54\phantom{\rule{thinmathspace}{0ex}}{s}^{1},{l}_{i}=0.52\phantom{\rule{thinmathspace}{0ex}}{s}^{1},\sigma =5.54\phantom{\rule{thinmathspace}{0ex}}{s}^{1}.$
The three alternative models produce freerunning oscillatory dynamics with similar frequency and amplitude as measured from worms swimming in fluids with viscosity $120\phantom{\rule{thinmathspace}{0ex}}mPa\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\phantom{\rule{thinmathspace}{0ex}}s$. Here, we modeled the process of muscle inhibition by multiplying the scaled active muscle moment, ${M}_{a}$, with a factor, $1Q\left(\mathrm{\Delta}t\right)$, as a function of the time interval $\mathrm{\Delta}t$ in a bellshaped 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}\left(t\right)$. Specifically, ${M}_{a}\left(t\right)$ is positive when ventral muscles contract and dorsal muscles relax, and negative for the other half of the cycle. Therefore, in our thresholdswitch model, specifically inhibiting dorsal or ventral or bothside muscles was computationally equivalent to conditionally modulating ${M}_{a}\left(t\right)$ with the bellshaped modulating function depending on the sign of ${M}_{a}\left(t\right)$.
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, $1Q\left(\mathrm{\Delta}t\right)$, which has a bellshaped profile (Figure 4—figure supplement 1):
where $r=0.3\phantom{\rule{thinmathspace}{0ex}}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\xb7{10}^{1/q}$ so that $Q{}_{\mathrm{\Delta}\mathrm{t}=0}>0.99$. Note that when modeling the dorsalsideonly 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 dorsalside illumination to show a relatively moderate response compared to ventralside illumination.
To simulate the muscle inhibition on our thresholdswitch model, we multiplied ${M}_{a}$ with $(1Q)$ 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 bellshaped 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 StuartLandau Oscillator, it is modified as:
For each individual model listed above, $\stackrel{~}{{M}_{i}}$ (subscript $i$ represents V, R, and S, respectively) is the generalized muscle moment which is to be multiplied by the bellshaped factor $(1Q)$ 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 A15A17 would be equivalent to their original form (Equations A11–A13) when inhibition is absent (in this case, $\stackrel{~}{{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 tworound fitting process. First, a subset of parameters was determined by fitting the model to observations of freemoving dynamics; the model could generate freemoving 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 twostep optimization procedure for individual models are provided as follows:
For the original thresholdswitch model, parameters ${\tau}_{u}$, ${M}_{0}$, and ${\tau}_{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 $\varphi $ with 100 evenly distributed samples. In this case, $\mathrm{\Delta}\varphi =2\pi /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 twostep optimization procedures for individual alternative models are summarized in Appendix 1—table 1.
Data availability
All data and software have been deposited to a Dryad repository (https://doi.org/10.5061/dryad.wwpzgmsk2).

Dryad Digital RepositoryData from: Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in Caenorhabditis elegans.https://doi.org/10.5061/dryad.wwpzgmsk2
References

The Kuramoto model: a simple paradigm for synchronization phenomenaReviews of Modern Physics 77:137–185.https://doi.org/10.1103/RevModPhys.77.137

Coordination and fine motor control depend on Drosophila TRPγNature Communications 6:1–13.https://doi.org/10.1038/ncomms8288

Peripheral feedback mechanisms acting on the central pattern generators for locomotion in fish and catCanadian Journal of Physiology and Pharmacology 59:713–726.https://doi.org/10.1139/y81108

The Caenorhabditis elegans unc49 locus encodes multiple subunits of a heteromultimeric GABA receptorThe Journal of Neuroscience 19:5348–5359.https://doi.org/10.1523/JNEUROSCI.191305348.1999

Sensory control of leg movement in the stick insect Carausius morosusBiological Cybernetics 25:61–72.https://doi.org/10.1007/BF00337264

Gait modulation in C. elegans: an integrated neuromechanical modelFrontiers in Computational Neuroscience 6:10.https://doi.org/10.3389/fncom.2012.00010

Neural control of Caenorhabditis elegans forward locomotion: the role of sensory feedbackBiological Cybernetics 98:339–351.https://doi.org/10.1007/s0042200802126

A consistent muscle activation strategy underlies crawling and swimming in Caenorhabditis elegansJournal of the Royal Society Interface 12:20140963.https://doi.org/10.1098/rsif.2014.0963

Sensory modification of leech swimming: interactions between ventral stretch receptors and swimrelated neuronsJournal of Comparative Physiology A 187:569–579.https://doi.org/10.1007/s003590100229

The neural circuit for touch sensitivity in Caenorhabditis elegansThe Journal of Neuroscience 5:956–964.https://doi.org/10.1523/JNEUROSCI.050400956.1985

BookThe Nematode Caenorhabditis elegansNew York: Cold Spring Harbor Laboratory Press.

Mechanisms underlying rhythmic locomotion: body–fluid interaction in undulatory swimmingJournal of Experimental Biology 214:561–574.https://doi.org/10.1242/jeb.048751

The neuronal correlate of locomotion in fishExperimental Brain Research 41:11–18.https://doi.org/10.1007/BF00236674

The behaviour of Nematodes, their activity, senses and responsesThe Journal of Applied Ecolog 8:956.https://doi.org/10.2307/2402696

Signatures of proprioceptive control in Caenorhabditis elegans locomotionPhilosophical Transactions of the Royal Society B: Biological Sciences 373:20180208.https://doi.org/10.1098/rstb.2018.0208

Biological rhythms considered as relaxation oscillationsActa Medica Scandinavica 103:76–88.https://doi.org/10.1111/j.09546820.1940.tb11082.x

BookMechanotransductionIn: Kaplan J, editors. C. elegans II. New York, United States: Cold Spring Harbor Laboratory Press. pp. 10–22.

Simulations of neuromuscular control in lamprey swimmingPhilosophical Transactions of the Royal Society of London. Series B: Biological Sciences 354:895–902.https://doi.org/10.1098/rstb.1999.0441

BookStatistical Analysis of Spherical DataCambridge University Press.https://doi.org/10.1017/CBO9780511623059

BookCentral pattern generators: sensory feedbackIn: Friesen WO, editors. Encyclopedia of Neuroscience. Elsevier. pp. 701–709.

Advances in Insect Physiology31–140, Pattern and control of walking in insects, Advances in Insect Physiology, Elsevier, 10.1016/S00652806(08)600399.

Studies in animal locomotion: I. the movement of fish with special reference to the eelJournal of Experimental Biology 10:88–104.https://doi.org/10.1242/jeb.10.1.88

The locomotion of nematodesJournal of Experimental Biology 41:135–154.https://doi.org/10.1242/jeb.41.1.135

The motor infrastructure: from ion channels to neuronal networksNature Reviews Neuroscience 4:573–586.https://doi.org/10.1038/nrn1137

Phaseresponse curves give the responses of neurons to transient inputsJournal of Neurophysiology 94:1623–1635.https://doi.org/10.1152/jn.00359.2004

Behavioral plasticity in C. elegans: paradigms, circuits, genesJournal of Neurobiology 54:203–223.https://doi.org/10.1002/neu.10168

From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegansPhilosophical Transactions of the Royal Society B: Biological Sciences 373:20170374.https://doi.org/10.1098/rstb.2017.0374

Neuromechanical Mechanisms of Gait Adaptation in C. elegans : Relative Roles of Neural and Mechanical CouplingSIAM Journal on Applied Dynamical Systems 20:1022–1052.https://doi.org/10.1137/20M1346122

Systems level circuit model of C. elegans undulatory locomotion: mathematical modeling and molecular geneticsJournal of Computational Neuroscience 24:253–276.https://doi.org/10.1007/s1082700700546

On the movement of slender bodies near plane boundaries at low Reynolds numberJournal of Fluid Mechanics 72:529–540.https://doi.org/10.1017/S0022112075003126

BookNeuronal Mechanisms for Generating Locomotor ActivityNew York, United States: New York Academy of Sciences.

Development and functional organization of spinal locomotor circuitsCurrent Opinion in Neurobiology 21:100–109.https://doi.org/10.1016/j.conb.2010.09.004

Rhythmic swimming activity in neurones of the isolated nerve cord of the leechJournal of Experimental Biology 65:643–668.https://doi.org/10.1242/jeb.65.3.643

Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotionPLOS Computational Biology 13:e1005303.https://doi.org/10.1371/journal.pcbi.1005303

Central pattern generators and the control of rhythmic movementsCurrent Biology 11:R986–R996.https://doi.org/10.1016/S09609822(01)005814

Principles of rhythmic motor pattern generationPhysiological Reviews 76:687–717.https://doi.org/10.1152/physrev.1996.76.3.687

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

Voltage oscillations in the barnacle giant muscle fiberBiophysical Journal 35:193–213.https://doi.org/10.1016/S00063495(81)847820

BookA PRC Description of How Inhibitory Feedback Promotes Oscillation stabilityIn: Bose A, editors. Phase Response Curves in Neuroscience. Springer. pp. 399–417.https://doi.org/10.1007/9781461407393_16

An active pulse transmission line simulating nerve axonProceedings of the IRE 50:2061–2070.https://doi.org/10.1109/JRPROC.1962.288235

A neuromechanical model of multiple network rhythmic pattern generators for forward locomotion in C. elegansFrontiers in Computational Neuroscience 15:10.https://doi.org/10.3389/fncom.2021.572339

BookGenerating the Walking Gait: Role of Sensory feedbackProgress in Brain ResearchAmsterdam, Netherlands: Elsevier.https://doi.org/10.1016/S00796123(03)430124

Flexible mechanisms: the diverse roles of biological springs in vertebrate movementJournal of Experimental Biology 214:353–361.https://doi.org/10.1242/jeb.038588

Kinetic properties of mechanically activated currents in spinal sensory neuronsThe Journal of Physiology 588:301–314.https://doi.org/10.1113/jphysiol.2009.182360

BookPhase Response Curves in NeuroscienceBerlin, Germany: Springer.https://doi.org/10.1007/9781461407393

Phaseresponse curves and synchronized neural networksPhilosophical Transactions of the Royal Society B: Biological Sciences 365:2407–2422.https://doi.org/10.1098/rstb.2009.0292

Scatter search and local NLP solvers: a multistart framework for global optimizationINFORMS Journal on Computing 19:328–340.https://doi.org/10.1287/ijoc.1060.0175

Lxxxviii.on “relaxationoscillations”The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2:978–992.https://doi.org/10.1080/14786442608564127

The motor circuitInternational Review of Neurobiology 69:125–167.https://doi.org/10.1016/S00747742(05)690058

The dynamics of nematode movementAnnual Review of Phytopathology 6:91–114.https://doi.org/10.1146/annurev.py.06.090168.000515

The structure of the nervous system of the nematode Caenorhabditis elegansPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 314:1–340.https://doi.org/10.1098/rstb.1986.0056

Sensory feedback can coordinate the swimming activity of the leechThe Journal of Neuroscience 19:4634–4643.https://doi.org/10.1523/JNEUROSCI.191104634.1999

Entrainment of leech swimming activity by the ventral stretch receptorJournal of Comparative Physiology A 190:939–949.https://doi.org/10.1007/s0035900405499
Decision letter

Ronald L CalabreseSenior Editor; Emory University, United States

Manuel ZimmerReviewing Editor; University of Vienna, Austria

Carter JohnsonReviewer; University of California, Berkeley, Berkeley, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This paper combines behaviour quantification, optogenetic perturbation, and modelling to study C. elegans locomotion. It will be of interest to neuroscientists studying lcomotion and gait adaptation. The quantitative agreement between the model and experimentsimportantly including the perturbation experimentssuggests forward locomotion in worms can be understood as being driven by a relaxation oscillator. This conclusion provides intuition for how worms move and should provide useful constraints on more detailed models that include neural anatomy and physiology.
Decision letter after peer review:
Thank you for submitting your article "Phase response analyses support a relaxation oscillator model of locomotor rhythm generation in Caenorhabditis elegans" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by Manuel Zimmer as a Reviewing Editor and Ronald Calabrese as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Carter Johnson (Reviewer #3).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. As you can see below, the reviewers provide an extensive list of detailed comments, however the overall result of the reviews and subsequent discussions is that we are very excited about your work, and we would like to support publication in eLife of a revised manuscripts that addresses all comments. The comments below are partially concerned with similar issues but to retain each reviewers' details, we provide you here with the full list. Please submit with your revision a pointbypoint response to the comments.
Essential revisions:Reviewer #1:
(R1.1) The authors set out to better understand the dynamics of C. elegans locomotion characterizing the rhythmic pattern of the worm's head during locomotion and by generating the phase response curve in that region. The authors then illustrate how well the relaxationoscillator phenomenological model, with only a handful of free parameters optimized to match the data, is consistent with the observations, including observations under different mechanical loads.
(R1.2) The authors provide a characterization of the curvature dynamics in the worm's head region (0.10.3 body coordinate). The reasoning for focusing on this area in particular is not discussed or justified in the manuscript. Also, no comment is made whether the characteristics observed for the limit cycle in that region is consistent across the rest of the regions. From the experimental description and data shown (i.e., Figure 2B), it would appear that the authors had access to a characterization of the limit cycle across any region of the body of the worm, from head to tail. To avoid the appearance of cherrypicking a region, authors are encouraged to include the full characterization of the limit cycle for all the data available to them. In addition to that, authors are encouraged to clarify explicitly whether the characteristics observed in one region (e.g., nonsinusoidal limit cycle) hold true across the rest of the worm's body. How sinusoidal or nonsinusoidal, nonellipsoidal, asymmetric is the limit cycle as a function of the body region of the worm, from head to tail? Crucially, if the characterization observed in the selected region is not equally pronounced across the rest of body, it seems crucial to discuss it in the interpretation of the results. This is particularly important given that throughout the manuscript, the authors interpret the results as applicable to the rhythmic generation across the whole body, not just that limited region (see first sentence of the Discussion, for example; and it is also relevant given that they discuss equally proprioceptive feedback in the head as well as in the rest of the body). These points apply to the characterization of the limit cycle given that the data is there but not shown; this is not the case for the phase response analysis, because understanding it for the full body would require further experimentation. However, this should be a point of discussion: is the expectation that the phase response curve through transient ontogenetic inhibition will look similar throughout the body? Is the same sawtoothshape with the sharp transitions in the same part of the phase expected? What is the implication of the expected results if they are consistent or not consistent with the ones shown for the region focused on here?
(R1.2) The model assumes proprioceptive feedback from the curvature of the same region that is bending. However, the undifferentiated processes that are used as inspiration for this assumption extend posteriorly in the case of SMD and anteriorly in the case of Bclass motor neurons. So the assumption that the curvature in one region can affect the bending in that same region should be discussed.
(R1.3) The model also assumes that the stretch receptors sense both the curvature and the rate of curvature. However, the authors provide no evidence of other C. elegans neurons where both of those can be sensed simultaneously. It seems feasible that a neuron might be sensitive to either the curvature or the rate of curvature, but it seems less feasible that the neuron would be sensitive to a linear combination of both. This seems to be a central assumption of the model, with little or no justification.
(R1.4) I commend the authors for comparing their selected model with other similar ones. I would highly encourage that the authors include a comparison of the models into the main manuscript. In particular, given the similarities between the results, I would highly encourage the authors to place the results of their chosen model next to the alternatives (it can be a figure like S10 but including the preferred model). As it is, with the alternatives in S10 and the rest of the figures from the best model in the main text, the comparison is significantly harder.
(R1.5) Finally, the discussion of previous relevant computational models of C. elegans locomotion is rather limited. There are at least four particularly relevant contributions that the authors do not discuss in relation to their own contribution:
(a) Johnson, Lewis, and Guy (2020) Neuromechanical Mechanisms of Gait Adaptation in C. elegans: Relative Roles of Neural and Mechanical Coupling.
(b) Izquierdo and Beer (2018) From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegans.
(c) Kunert, Proctor, Brunton, Kutz (2017) Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotion
(d) Denham, Ranner, Cohen (2018) Signatures of proprioceptive control in Caenorhabditis elegans locomotion.
Reviewer #2:
(R2.1) The paper's conclusions are supported by the data and the experiments and model results are clearly presented. In some cases there could be more analysis of uncertainty. For example, how consistent is the phase portrait? This could be illustrated by including some individual traces and/or by estimating a confidence interval (using resampling if necessary). A different kind of uncertainty that should be explained is that arising from fitting the model. Was the optimisation sensitive to the choice of starting parameters? Are there several local minima that explain the data reasonably well or did the optimisation always converge to the same parameter values?
(R2.2) The authors have focused on forward locomotion bouts and have also eliminated cycles with a period that deviated more than 20% from the average across worms within a session. How many of the bouts were eliminated with this threshold? Are any of the conclusions of the paper sensitive to this threshold?
(R2.3) The authors have convincingly ruled out the van der Pol and StuartLandau oscillators but the Rayleigh oscillator fits the data quite well with three parameters. Some quantification of the relative fit quality of the thresholdswitch model and the Rayleigh might help guide future work.
(R2.4) Finally, the discussion does a good job of putting the results of the phenomenological thresholdswitch model in the context of previously published detailed models but it might be strengthened with more quantitative comparisons. For example, do the observed muscle switching time scales or the nature of the proprioceptive threshold rule out or usefully constrain any of the detailed models?
(R2.5) The schematic figures (1 and 4) would be improved by leaving out the abbreviations. At a minimum, include (BWM) after 'body wall muscles' in the caption of figure 1.
(R2.6) Would a different colour map make Figure S2 clearer? Jet makes some regions stand out more than is warranted. For example, the purple patch in the lower right corner and red patch in the upper left corder are visually very salient but are probably just normal fluctuations within the noise.
(R2.7) In Figure 3, grey is used to indicate the control data in B but is used for individual trials in subsequent panels. Use a different colour for controls in B? Also add the grey curves to the legend next to D rather than only describing them in the caption.
(R2.8) Line 299, state the gene name of the GABA receptor.
(R2.9) In the Appendix there are a few cases of missing articles (for example line 812, "worm's head" should be "the worm's head" and on line 886 "animal's" should be "the animal's").
(R2.10) Starting on line 817, D^C is used to denote the distance from the origin over the normal cycle. To avoid confusion with the scaling constant c, perhaps choose a different symbol for the scaling constant or to indicate the normal cycle?Reviewer #3:
(R3.1) The main suggestion I have is to change the framing of the section "The Analysis of Alternative Models Supports a Relaxation Oscillation Mechanism'. In the first paragraph, the authors state that the purpose of this section is to "[ask] whether other models could also explain our findings." However, the section title and final paragraph suggest that the purpose of the section is to support the idea that the underlying oscillatory mechanism is a relaxationoscillation. If my understanding of the authors' purpose here is correct, then I don't think the analysis of these three specific phenomenological models provides evidence of this. Specifically, the last paragraph of this section (Lines 586590), suggests that because the nonrelaxation oscillator (StuartLandau) was not able to capture your results (nonsinusoidal limit cycle and sawtooth PRCs), the underlying mechanism is a relaxationoscillation. I don't think these models support that claim, because a different nonrelaxation oscillator tuned precisely may be able to capture these results. I think a much more thorough analysis of general oscillator dynamics would be needed to make such a claim. To be clear, I think the paper is strong enough without this analysis, I would just be more cautious about the claims here.
(R3.2) The strength of the author's original computational model is that it is based on specific mechanisms thought to underlie C. elegans neurolocomotion (specifically motor neurons, muscles, and proprioception). The authors discuss how this model functions as a relaxationoscillator in the discussion, but I think it would be stronger in this section.
(R3.3) Intro or discussion – In addition to Cohen's group's models, two recent modeling papers from other modeling groups have investigated the C. elegans neurolocomotion system explicitly as a chain of coupled neuromechanical oscillators (Olivares et al., 2021 doi:10.3389/fncom.2021.572339 and Johnson et al., 2021 doi:10.1137/20M1346122). In particular, Johnson et al., 2021 show how phaseresponse properties influence coordination of the neuromechanical oscillators. Describing how your computational model compares with or supports these other recent models would help give context to the modeling contribution of this work.
(R3.14) Lines 311313: The authors mention that sawtooth PRCS "may reflect a phaseresetting property of an oscillator with respect to a perturbation". I think more detail would be helpful here. To my knowledge, this is referring to the idea that if an oscillator gets a really strong perturbation, the phaseresponse curve turns into a "phaseresetting curve", where the large perturbation essentially resets the phase of the oscillator, and then it just marches forward linearly through time. Are you suggesting that this sawtooth shape is indicative of your perturbation being large? This would be an important detail to include.
(R3.5) Section "Analysis of Alternative Models Supports a Relaxation Oscillation Mechanism":
The strength of the author's original computational model is that it is based on specific mechanisms thought to underlie C. elegans neurolocomotion (specifically motor neurons, muscles, and proprioception). To support the idea that the underlying mechanism is a relaxationoscillation, I would suggest that the authors instead analyze their own computational model and show or point out how it functions structurally as a relaxationoscillation. This is mentioned briefly in the discussion, but I think it would be stronger in this section. Furthermore, I think the authors could give the reader a clearer reason why exactly showing that the mechanism is a relaxationoscillation is important. To my understanding, the open question in the C. elegans locomotionmodeling literature is whether the oscillations are fundamentally neurallydriven (like a CPG or HCO) or reflexdriven (like the proprioceptive mechanism here). Perhaps evidence for the relaxationoscillation mechanism could be made by considering an alternative mechanistic model with neurallydriven oscillations and whether it can capture the new phase response data.
(R3.6) The authors have not given adequate description of the model fitting process in the main text. In line 548, you only mention that appropriate parameters were chosen, and reference the supplement. I don't think too much detail is needed here, but explaining that the models were fit to match both the limit cycles and PRCs would be helpful here. My first impression of this section was that the limit cycles were matched and the PRCs were emergent. After looking at the supplementary details I saw that you explicitly attempted to match both results. This would be helpful to put up front.
(R3.7) Lines 6164: Sentence grammar/structure. The "but also …" is missing a "not only" earlier in the sentence (e.g. "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; 13).").
(R3.8) Lines 152153 and Figure 2D caption: I think it would help to clarify in either or both of these places that the boxed portion of the nematode bodies corresponds to the curvature region (0.10.3 body coordinates).
(R3.9) Lines 940951: The mechanics here depend on the assumption of a traveling wave of fixed wavelength λ down the body. I would mention explicitly that both the wavelength λ and normal drag coefficient CN will be varied with fluid viscosity.
https://doi.org/10.7554/eLife.69905.sa1Author response
Essential revisions:
Reviewer #1:
(R1.1) The authors set out to better understand the dynamics of C. elegans locomotion characterizing the rhythmic pattern of the worm's head during locomotion and by generating the phase response curve in that region. The authors then illustrate how well the relaxationoscillator phenomenological model, with only a handful of free parameters optimized to match the data, is consistent with the observations, including observations under different mechanical loads.
(R1.2) The authors provide a characterization of the curvature dynamics in the worm's head region (0.10.3 body coordinate). The reasoning for focusing on this area in particular is not discussed or justified in the manuscript. Also, no comment is made whether the characteristics observed for the limit cycle in that region is consistent across the rest of the regions. From the experimental description and data shown (i.e., Figure 2B), it would appear that the authors had access to a characterization of the limit cycle across any region of the body of the worm, from head to tail. To avoid the appearance of cherrypicking a region, authors are encouraged to include the full characterization of the limit cycle for all the data available to them. In addition to that, authors are encouraged to clarify explicitly whether the characteristics observed in one region (e.g., nonsinusoidal limit cycle) hold true across the rest of the worm's body. How sinusoidal or nonsinusoidal, nonellipsoidal, asymmetric is the limit cycle as a function of the body region of the worm, from head to tail? Crucially, if the characterization observed in the selected region is not equally pronounced across the rest of body, it seems crucial to discuss it in the interpretation of the results. This is particularly important given that throughout the manuscript, the authors interpret the results as applicable to the rhythmic generation across the whole body, not just that limited region (see first sentence of the Discussion, for example; and it is also relevant given that they discuss equally proprioceptive feedback in the head as well as in the rest of the body).
As suggested, we have calculated the limit cycles of bending dynamics as a function of body coordinate from the head to the tail and included the results in the manuscript. We also provide a justification for focusing on the head region for most of our analyses (see Lines 159162 and Figure 2—figure supplement 1).
These points apply to the characterization of the limit cycle given that the data is there but not shown; this is not the case for the phase response analysis, because understanding it for the full body would require further experimentation. However, this should be a point of discussion: is the expectation that the phase response curve through transient ontogenetic inhibition will look similar throughout the body? Is the same sawtoothshape with the sharp transitions in the same part of the phase expected? What is the implication of the expected results if they are consistent or not consistent with the ones shown for the region focused on here?
We conducted further experiments to optogenetically inhibit body wall muscles in other regions of a worm besides the head region. The corresponding PRCs are now shown and discussed in the manuscript (see Lines 265269, and Lines 952978 and Figure 3—figure supplement 5 in Appendix).
(R1.2) The model assumes proprioceptive feedback from the curvature of the same region that is bending. However, the undifferentiated processes that are used as inspiration for this assumption extend posteriorly in the case of SMD and anteriorly in the case of Bclass motor neurons. So the assumption that the curvature in one region can affect the bending in that same region should be discussed.
The thresholdswitch model as well as the other models are developed based on experimental observations of the head region, defined as the body coordinate interval between 0.1 and 0.3, which is about 200 μm long. The asynaptic processes in the posterior portion of the SMD neurons are ~40 μm long (c.f. the graphical abstract for Reid et al., Neuroscience 2015, PMID: 26480814) and have been covered by the region of our analyses. The asynaptic processes in the anterior portion of the nearhead Btype motor neurons are 100~200 μm long (see Figure 1 in Chen et al., PNAS 2006, PMID: 16537428) where most of the processes have also been covered by the region of our analyses. Therefore, the processes of the candidate proprioceptors (SMD and anterior Btype motor neurons) could affect the bending in the analyzed region.
(R1.3) The model also assumes that the stretch receptors sense both the curvature and the rate of curvature. However, the authors provide no evidence of other C. elegans neurons where both of those can be sensed simultaneously. It seems feasible that a neuron might be sensitive to either the curvature or the rate of curvature, but it seems less feasible that the neuron would be sensitive to a linear combination of both. This seems to be a central assumption of the model, with little or no justification.
First, it does not seem at all implausible to us that a C. elegans neuron might be sensitive to a combination of curvature and rate of change of curvature. For example, in rat dorsal root ganglia, rapidlyadapting (RA) and intermediateadapting (IA) mechanosensors encode a combination of stimulus size and velocity. See for example Figure 1 in Rugiero et al., J Physiol 2010 (PMID: 19948656).
Second, even if we assume a single neuron is not sensitive to a combination of curvature and rate of change of curvature, a collection of neurons involved in motor rhythm generation may be. The stretch receptors in our model might be considered as a collective behavior of curvaturesensing neurons and curvaturechangesensing
neurons.
We have revised this discussion to clarify these points (see Lines 714729).
(R1.4) I commend the authors for comparing their selected model with other similar ones. I would highly encourage that the authors include a comparison of the models into the main manuscript. In particular, given the similarities between the results, I would highly encourage the authors to place the results of their chosen model next to the alternatives (it can be a figure like S10 but including the preferred model). As it is, with the alternatives in S10 and the rest of the figures from the best model in the main text, the comparison is significantly harder.
As suggested, we have added our primary model to the comparison between all models (see Figure 5—figure supplement 1 and Figure 6—figure supplement 2).
(R1.5) Finally, the discussion of previous relevant computational models of C. elegans locomotion is rather limited. There are at least four particularly relevant contributions that the authors do not discuss in relation to their own contribution:
(a) Johnson, Lewis, and Guy (2020) Neuromechanical Mechanisms of Gait Adaptation in C. elegans: Relative Roles of Neural and Mechanical Coupling.
(b) Izquierdo and Beer (2018) From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegans.
(c) Kunert, Proctor, Brunton, Kutz (2017) Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotion
(d) Denham, Ranner, Cohen (2018) Signatures of proprioceptive control in Caenorhabditis elegans locomotion.
We thank the reviewer for these suggestions. These models have been cited and discussed in the revised manuscript (see Lines 117119, 639640, 646648, and 672674).
Reviewer #2:
(R2.1) The paper's conclusions are supported by the data and the experiments and model results are clearly presented. In some cases there could be more analysis of uncertainty. For example, how consistent is the phase portrait? This could be illustrated by including some individual traces and/or by estimating a confidence interval (using resampling if necessary).
As suggested, we have added 10 randomly selected individual traces to the average
trace shown in Figure 2D.
A different kind of uncertainty that should be explained is that arising from fitting the model. Was the optimisation sensitive to the choice of starting parameters? Are there several local minima that explain the data reasonably well or did the optimisation always converge to the same parameter values?
Our fitting procedure for optimizing model parameters was not sensitive to the initial choices. For each individual model developed in this manuscript, the parameters were determined via a tworound fitting procedure as discussed in Results (Lines 367373) and Appendix (Lines 12171248). In the first round, a model’s parameters were optimized to match the experimental free locomotion, using cost functions shown in Equation S7 (for thresholdswitch model) and Table S1 (for alternative models). During this round, the parameters converged to a unique set of values since there is only one local minimum within the range of parameters (a large range of parameters that could generate oscillation). In the second round, the remaining undetermined parameters were optimized to match the experimental PRC curve using a global minimum search strategy within a range of parameters. Since it is a global search algorithm, the optimization is independent of the choice of initial parameters.
(R2.2) The authors have focused on forward locomotion bouts and have also eliminated cycles with a period that deviated more than 20% from the average across worms within a session. How many of the bouts were eliminated with this threshold?
Under the threshold 20%, 857 out of 1910 trials (~45%) that deviated from the average period were eliminated.
Are any of the conclusions of the paper sensitive to this threshold?
None of the conclusions of the paper are sensitive to the threshold. We measured the PRC and limit cycle under thresholds up to 50% (343 out of 1910 bouts (~18%) were dropped), and down to 10% (1270 out of 1910 bouts (67%) were dropped). We find that the resulting curves are largely independent of the choice of threshold (see Author response image 1).
(R2.3) The authors have convincingly ruled out the van der Pol and StuartLandau oscillators but the Rayleigh oscillator fits the data quite well with three parameters. Some quantification of the relative fit quality of the thresholdswitch model and the Rayleigh might help guide future work.
As now shown in Figure 5—figure supplement 1, the thresholdswitch model fits the experimental PRC better than Rayleigh model by a considerable margin in terms of mean square error (thresholdswitch: 0.12; vdP: 0.21; Rayleigh: 0.37).
(R2.4) Finally, the discussion does a good job of putting the results of the phenomenological thresholdswitch model in the context of previously published detailed models but it might be strengthened with more quantitative comparisons. For example, do the observed muscle switching time scales or the nature of the proprioceptive threshold rule out or usefully constrain any of the detailed models?
As suggested, we have added an extensive discussion to acknowledge the potential usefulness of parameters measurements for constraining detailed models (see Lines 660674).
(R2.5) The schematic figures (1 and 4) would be improved by leaving out the abbreviations. At a minimum, include (BWM) after 'body wall muscles' in the caption of figure 1.
We have added the full names for all the abbreviations shown in the Figures (1 and 4) and the corresponding captions.
(R2.6) Would a different colour map make Figure S2 clearer? Jet makes some regions stand out more than is warranted. For example, the purple patch in the lower right corner and red patch in the upper left corder are visually very salient but are probably just normal fluctuations within the noise.
We have changed this figure to have a cyclic color map with a constant lightness which avoids visual emphasis on certain parts of the fluctuations (see Figure 3— figure supplement 2).
(R2.7) In Figure 3, grey is used to indicate the control data in B but is used for individual trials in subsequent panels. Use a different colour for controls in B? Also add the grey curves to the legend next to D rather than only describing them in the caption.
As suggested, we have revised Figure 3 and its caption.
(R2.8) Line 299, state the gene name of the GABA receptor.
We have added the gene information in the main text (Lines 307309) and the figure captions (Figure 3—figure supplement 8). Also, strain names that previously appeared in figure captions have been moved to Methods.
(R2.9) In the Appendix there are a few cases of missing articles (for example line 812, "worm's head" should be "the worm's head" and on line 886 "animal's" should be "the animal's")
We have corrected the multiple missing articles pointed out by the reviewer.
(R2.10) Starting on line 817, D^C is used to denote the distance from the origin over the normal cycle. To avoid confusion with the scaling constant c, perhaps choose a different symbol for the scaling constant or to indicate the normal cycle?
We have used a different symbol for the scaling constant (see Figure 3—figure supplement 2 and Lines 838867).
Reviewer #3:
(R3.1) The main suggestion I have is to change the framing of the section "The Analysis of Alternative Models Supports a Relaxation Oscillation Mechanism'. In the first paragraph, the authors state that the purpose of this section is to "[ask] whether other models could also explain our findings." However, the section title and final paragraph suggest that the purpose of the section is to support the idea that the underlying oscillatory mechanism is a relaxationoscillation. If my understanding of the authors' purpose here is correct, then I don't think the analysis of these three specific phenomenological models provides evidence of this. Specifically, the last paragraph of this section (Lines 586590), suggests that because the nonrelaxation oscillator (StuartLandau) was not able to capture your results (nonsinusoidal limit cycle and sawtooth PRCs), the underlying mechanism is a relaxationoscillation. I don't think these models support that claim, because a different nonrelaxation oscillator tuned precisely may be able to capture these results. I think a much more thorough analysis of general oscillator dynamics would be needed to make such a claim. To be clear, I think the paper is strong enough without this analysis, I would just be more cautious about the claims here.
We agree that our investigation of a few alternative models may not lead to the conclusion claimed by the original section title. We have rephrased this section and changed its title appropriately (see Lines 551608).
(R3.2) The strength of the author's original computational model is that it is based on specific mechanisms thought to underlie C. elegans neurolocomotion (specifically motor neurons, muscles, and proprioception). The authors discuss how this model functions as a relaxationoscillator in the discussion, but I think it would be stronger in this section.
See response to R3.5.
(R3.3) Intro or discussion – In addition to Cohen's group's models, two recent modeling papers from other modeling groups have investigated the C. elegans neurolocomotion system explicitly as a chain of coupled neuromechanical oscillators (Olivares et al., 2021 doi:10.3389/fncom.2021.572339 and Johnson et al., 2021 doi:10.1137/20M1346122). In particular, Johnson et al., 2021 show how phaseresponse properties influence coordination of the neuromechanical oscillators. Describing how your computational model compares with or supports these other recent models would help give context to the modeling contribution of this work.
As suggested, we have discussed and cited these two papers in the Introduction (Lines 117119) and Discussion (Lines 639641, 646648, and 672674).
(R3.14) Lines 311313: The authors mention that sawtooth PRCS "may reflect a phaseresetting property of an oscillator with respect to a perturbation". I think more detail would be helpful here. To my knowledge, this is referring to the idea that if an oscillator gets a really strong perturbation, the phaseresponse curve turns into a "phaseresetting curve", where the large perturbation essentially resets the phase of the oscillator, and then it just marches forward linearly through time. Are you suggesting that this sawtooth shape is indicative of your perturbation being large? This would be an important detail to include.
We did not intend the term phase resetting to be different from phase response in our manuscript. We apologize for any confusion and have revised this passage.
To the main point of the question, we did also perform the PRC experiments with weaker stimuli (0.055 s muscle inhibition, Figure 3—figure supplement 4), which resulted in a similar sawtoothshaped PRC.
(R3.5) Section "Analysis of Alternative Models Supports a Relaxation Oscillation Mechanism":
The strength of the author's original computational model is that it is based on specific mechanisms thought to underlie C. elegans neurolocomotion (specifically motor neurons, muscles, and proprioception). To support the idea that the underlying mechanism is a relaxationoscillation, I would suggest that the authors instead analyze their own computational model and show or point out how it functions structurally as a relaxationoscillation. This is mentioned briefly in the discussion, but I think it would be stronger in this section.
In Discussion (Lines 675680 and 694702), we have explained how our model satisfies the properties of a relaxation oscillator. In the caption of Figure 4A and Appendix (Lines 10181030), each site respectively contains a stepbystep explanation of how our model is constructed as a relaxation oscillator.
Furthermore, I think the authors could give the reader a clearer reason why exactly showing that the mechanism is a relaxationoscillation is important. To my understanding, the open question in the C. elegans locomotionmodeling literature is whether the oscillations are fundamentally neurallydriven (like a CPG or HCO) or reflexdriven (like the proprioceptive mechanism here). Perhaps evidence for the relaxationoscillation mechanism could be made by considering an alternative mechanistic model with neurallydriven oscillations and whether it can capture the new phase response data.
As highlighted in Introduction/Discussion (Lines 120125, 655659, and 725729) and described in Results/Appendix (Figure 4A, Lines 349373 and 9851041), our model is designed to be a phenomenological model that does not specify the neural and synaptic details of the motor circuit but rather capture the essential mechanisms underlying motor rhythm generation. We have compared our model with several alternatives based on other types of oscillators. We feel that consideration of alternative models with neurallydriven oscillations would be interesting but beyond the scope of the work.
In our view, the importance of showing that C. elegans locomotion operates by a relaxationoscillation mechanism is reflected by demonstrating how and why such a simple model with a few parameters could capture several different experimental observations (freelocomotion, phase responses, and gait adaptation to external load), and in the potential to obtain an integrative understanding of how an organism
generates behavior.
(R3.6) The authors have not given adequate description of the model fitting process in the main text. In line 548, you only mention that appropriate parameters were chosen, and reference the supplement. I don't think too much detail is needed here, but explaining that the models were fit to match both the limit cycles and PRCs would be helpful here. My first impression of this section was that the limit cycles were matched and the PRCs were emergent. After looking at the supplementary details I saw that you explicitly attempted to match both results. This would be helpful to put up front.
As suggested, we have moved the tworound fitting procedure to the beginning of the corresponding Result section to better clarify this process (see Lines 369373). More details can also be found in the Appendix (Lines 12171248)
(R3.7) Lines 6164: Sentence grammar/structure. The "but also …" is missing a "not only" earlier in the sentence (e.g. "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; 13).").
Thanks for pointing out this grammatical error. We have corrected this.
(R3.8) Lines 152153 and Figure 2D caption: I think it would help to clarify in either or both of these places that the boxed portion of the nematode bodies corresponds to the curvature region (0.10.3 body coordinates).
As suggested, we have clarified the contents in both places that the boxed region indicates the anterior region of the worm body (0.10.3 body coordinates).
(R3.9) Lines 940951: The mechanics here depend on the assumption of a traveling wave of fixed wavelength λ down the body. I would mention explicitly that both the wavelength λ and normal drag coefficient CN will be varied with fluid viscosity.
As suggested, we have modified this subsection to acknowledge that the derivation of the mechanics depends on the assumption of a travelling wave with a fixed wavelength (see Lines 10041017). The dependence of the wavelength and drag coefficient on the viscosity has been noted after the derivation.
https://doi.org/10.7554/eLife.69905.sa2Article and author information
Author details
Funding
National Institutes of Health (1R01NS084835)
 Hongfei Ji
 Anthony D Fouad
 Christopher FangYen
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 CFY. were supported by the National Institutes of Health (1R01NS084835). ST was supported by an Abraham Noordergraaf Research Fellowship and a Littlejohn Fellowship.
Senior Editor
 Ronald L Calabrese, Emory University, United States
Reviewing Editor
 Manuel Zimmer, University of Vienna, Austria
Reviewer
 Carter Johnson, University of California, Berkeley, Berkeley, United States
Publication history
 Preprint posted: June 23, 2020 (view preprint)
 Received: April 29, 2021
 Accepted: September 24, 2021
 Accepted Manuscript published: September 27, 2021 (version 1)
 Version of Record published: November 1, 2021 (version 2)
Copyright
© 2021, Ji et al.
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.
Metrics

 1,109
 Page views

 140
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Neuroscience
Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for highlevel tasks such as action recognition. To test this theory, we pursued a taskdriven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and modelbased muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the endeffector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the actionrecognitiontrained models and not the trajectorydecodingtrained models. This suggests that proprioceptive encoding is additionally associated with higherlevel functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.

 Computational and Systems Biology
Correlation between objects is prone to occur coincidentally, and exploring correlation or association in most situations does not answer scientific questions rich in causality. Causal discovery (also called causal inference) infers causal interactions between objects from observational data. Reported causal discovery methods and singlecell datasets make applying causal discovery to single cells a promising direction. However, evaluating and choosing causal discovery methods and developing and performing proper workflow remain challenges. We report the workflow and platform CausalCell (http://www.gaemons.net/causalcell/causalDiscovery/) for performing singlecell causal discovery. The workflow/platform is developed upon benchmarking four kinds of causal discovery methods and is examined by analyzing multiple singlecell RNAsequencing (scRNAseq) datasets. Our results suggest that different situations need different methods and the constraintbased PC algorithm with kernelbased conditional independence tests work best in most situations. Related issues are discussed and tips for best practices are given. Inferred causal interactions in single cells provide valuable clues for investigating molecular interactions and gene regulations, identifying critical diagnostic and therapeutic targets, and designing experimental and clinical interventions.