Neuromechanical Mechanisms of Gait Adaptation in C. elegans: Relative Roles of Neural and Mechanical Coupling

To simulate the model described in Section 2.1, the body is discretized into six modules in correspondence with the six neuromuscular modules, so that there are six discrete body segment curvatures. The 4th-order difference operator $D_{4}$ is used to approximate the 4th spatial derivative with zero-force, zero-torque boundary conditions:

Discretizing equations 3-5 and using 8-9 yields a linear differential equation for the vector of body segment curvatures $\underline{\mathbf{\kappa}}$:

where $I_{6}$ is the $6\times 6$ identity matrix.

In this discretization, the neural and muscle activity dynamics of all the modules are given by

where each vector entry (e.g., $A_{V,j}$) is the corresponding activity of the $j^{\text{th}}$ module. In equations 20 and 21, $W_{p}$ is the nonlocal proprioceptive connectivity matrix (equation 22), which comes from discretizing equation 13, and $W_{g}$ is the gap-junction connectivity matrix (equation 23), which comes from discretizing equation 15:

Some parameters in the model are well-constrained by experimental data, while others are not. Quantities that are directly measurable include the body length $L=1$ mm, average body radius $R=40$ $\mu$m, cuticle width $r_{c}=0.5$ $\mu$m, and wavelength $\lambda/L$ and frequency $f$ in fluids of various viscosities $\mu_{f}$. The timescales in the system are less certain. The range 50-200 ms is used for the muscle activation timescale $\tau_{m}$, which is the range of measurements of peak muscle force generation time in Milligan et al. (1997) [27]. As with previous models [5, 12, 19], the neural activity is chosen to be the fastest process in the model, but while Boyle et al. [5] considered the B-neurons as instantaneous switches, here the neural activity timescale is set at $\tau_{n}=10$ ms.

The internal viscosity $\mu_{b}$ and Young’s modulus $E$ have been estimated across several orders of magnitude [2, 13, 40], so caution is exercised in using one set of parameters from one source over another. Of more importance in the model is the mechanical timescale

which is the timescale of relaxation in an inviscid fluid. In equation 24, $k_{b}$ is the bending modulus, which is derived from the Young’s modulus $E$ and the geometry of the cuticle in Appendix B following previous modeling procedures [10, 40]. In Appendix B, we also compare the range of mechanical parameters used here with previous models. Given the range of mechanical parameters reported in the literature, the mechanical timescale could be as small as $\tau_{b}=1$ ms or as large as $\tau_{b}=5$ s. The role of this timescale is explored in Section 3.2.

The electrophysiological details of the internal neural circuit are largely unknown, thus all the feedback and coupling strengths $c_{p},c_{m},\varepsilon_{p},\varepsilon_{g}$, the parameters of the nonlinear functions $F(V)$ and $\sigma(A)$ are not well constrained. The feedback strengths $c_{m}=10$, $c_{p}=1$ and parameters of the nonlinear functions $F(V)$ ($a=1,I=0$) and $\sigma(A)$ ($c_{s}=1,a_{0}=2$) were chosen so that the neuromechanical oscillator robustly gives the correct frequency ($\sim 1.76Hz$) in a low-viscosity environment. The values for the coupling parameters $\varepsilon_{p}$ and $\varepsilon_{g}$, on the other hand, are explored in the next section.

We fit the model to match the wavelength and frequency in water, and then ran it in different fluid environments. Our model captures the quantitative effect of external fluid viscosity on the body wavelength seen in experiments and previous models. Figure 2 shows an example of the wavelength trend of the model for fixed body parameters $\tau_{b}=500$ ms, $\tau_{m}=50$ms (the wavelengths were computed from the model output as described in Appendix A). Figure 2 also shows that the model wavelengths are in close quantitative agreement with the experimentally-measured wavelengths of Fang-Yen et al. [13]. In water ($\mu_{f}=1$ mPa$\cdot$s) the wavelength is roughly 1.5 bodylengths, and increasing the fluid viscosity $\mu_{f}$ smoothly reduces the wavelength down to roughly 0.75 bodylengths in the most viscous case ($\mu_{f}=2.8\times 10^{4}$ mPa$\cdot$s). This wavelength trend is similar to what has been observed in other experiments [3, 40], and in Section 3.2, we show that our model captures this trend robustly over a wide range of parameters.

The undulation frequency also changes in response to changes in the fluid viscosity [3, 13, 40]. In Fang-Yen et al. [13], the frequency decreases from 1.7 Hz to 0.30 Hz as fluid viscosity increases from 1 mPa s to $2.8\times 10^{4}$ mPa s. Our model also exhibits a decrease in frequency as fluid viscosity $\mu_{f}$ increases, but not of the same magnitude (1.7 Hz - 1.6 Hz). Discussion of this discrepancy is given in Section 5.

We performed a parameter study to show that the model robustly captures gait adaptation as the fluid viscosity $\mu_{f}$ is varied. The mechanical parameters $\tau_{b}$, $\mu_{b}$, the proprioceptive coupling strength $\varepsilon_{p}$, and the muscle timescale $\tau_{m}$ were varied, while the other parameters of the model were held fixed, including the gap-junctional coupling strength $\varepsilon_{g}=0.0134$. (For more extensive parameter explorations, see [22].) For some parameter regimes, the body deformations were traveling waves for all fluid viscosities $\mu_{f}$, but this was not the case for other parameter regimes. Figure 3 shows kymographs of the body curvature that demonstrate two typical cases exhibited by the model.

The model parameters $\tau_{b}$, $\mu_{b}$, $\varepsilon_{p},$ and $\tau_{m}$ were systematically varied to characterize the model behavior. For a given body timescale $\tau_{b}$ and body viscosity $\mu_{b}$, the muscle activity timescale $\tau_{m}$ was selected in the range 50-250 ms to match the undulation frequency (1.7 Hz) in water ($\mu_{f}=1$ mPa s) within 1%. Next, the proprioceptive coupling strength $\varepsilon_{p}$ was selected to match the wavelength (1.5 bodylengths) in water within 1%. We are able to separate these effects and perform these one-parameter searches due to the weakly coupled nature of our model. The model was then run in different fluid viscosity $\mu_{f}$ environments and the emergent coordination trend is reported in Figure 4. The model behavior was classified exclusively as either: (1) not a traveling wave for all fluid environments, (2) incorrect wavelength trend, (3) qualitatively correct wavelength trend, or (4) incorrect frequency in water.

There is no traveling wave (red triangles) if, for any viscosity $\mu_{f}$, the difference between the minimum and maximum pairwise-phase difference is greater than or equal to 0.5, because this indicates that there is no consistent directionality to the phase differences in the body. A range of observed wavelength trends in various parameter regimes (the boxed markers in Figure 4) are illustrated in Figure 5. Figure 5(a) and (b) show examples of the qualitatively correct wavelength trend (blue circles), while (c) shows the the incorrect trend, which was only obtained at a single parameter combination. The wavelength trend is incorrect because the wavelengths increased dramatically as the fluid viscosity increased, as opposed to generally decreasing.

A few key observations can be made from Figure 4. First, if the mechanical timescale $\tau_{b}$ is too large, then the frequency in water cannot be obtained (see the black squares in Figure 4). Second, if the mechanical timescale $\tau_{b}$ is too small, then there will not be a traveling wave for all fluid viscosities $\mu_{f}$. This suggests that while the body stiffness $k_{b}$ and body viscosity $\mu_{b}$ have been estimated across several orders of magnitude in various experiments and models, the effective mechanical body timescale $\tau_{b}=\mu_{b}/k_{b}$ lies within the relatively narrow range $0.07-1$ s.

In order to match the frequency, the muscle timescale $\tau_{m}$ must be inversely related to $\tau_{b}$. When the body timescale $\tau_{b}$ is increased, the muscle timescale $\tau_{m}$ must decrease to compensate. The frequency in water cannot be obtained for $\tau_{b}$ too large since it would require decreasing the muscle activity timescale $\tau_{m}$ below physiological limits. Similarly, when the body timescale $\tau_{b}$ is decreased, the muscle timescale $\tau_{m}$ must be increased to compensate for the frequency. For $\tau_{b}$ too small, there is not a traveling wave for all fluid viscosities $\mu_{f}$; this occurs soon after $\tau_{b}<\tau_{m}$. This suggests that the relative ordering of the timescales $\tau_{b},\tau_{m},\tau_{n}$ is key to the coordination. Generally, the mechanical timescale $\tau_{b}$ must be the largest, the muscle activity timescale $\tau_{m}$ intermediate, and the neural timescale $\tau_{n}$ the shortest. The mechanism by which this timescale ordering affects coordination is explained in Section 4.3.

Remarkably, whenever there is a traveling wave in this systematic parameter search, it almost always has the qualitatively correct wavelength trend. This wavelength trend is consistent with gait adaptation across several orders of magnitude of the mechanical parameters.

To determine the role of gap-junctions in our neuromechanical model, we investigate how changing the gap-junctional coupling strength affects the wavelength trend and overall coordination. Other models have achieved accurate wavelength adaptation without gap-junctional coupling, albeit with proprioceptive ranges longer than nearest-neighbor [5, 12], so we assess the degree to which gap-junctional coupling is necessary in our model. We find that in our model, sufficiently strong gap-junctional coupling is necessary to obtain the long wavelengths in low-viscosity environments (e.g. water).

We vary the gap-junctional coupling strength in our model over several orders of magnitude ($\varepsilon_{g}=0$ and $\varepsilon_{g}\in[10^{-4},10^{-1}]$). We fix the parameters $\tau_{b}=0.5$ s, $\tau_{m}=0.1$ s, and $\mu_{b}=1.3\times 10^{-7}$ N$(mm)^{2}$s and find the proprioceptive coupling strength $\varepsilon_{p}$ to match the wavelength in water (approximately 1.5 bodylengths). Figure 6(a) shows the resulting wavelength trend of the model for different gap-junctional coupling strengths $\varepsilon_{g}$. For sufficiently strong gap-junctional coupling strengths ($\varepsilon_{g}\geq 0.005$), the wavelength in water was matched and the full range of wavelength adaptation was found (circle-markers in Figure 6(a)). For weaker gap-junctional coupling strengths ($\varepsilon_{g}<0.005$), the wavelength in water was not matched (plus-markers, Figure 6(a)). We further show that the model with weak gap-junctional coupling strength ($\varepsilon_{g}<0.005$) cannot match the wavelength in water for any proprioceptive coupling strength. Figure 6(b) shows the model wavelengths in water as a function of proprioceptive coupling strength $\varepsilon_{p}$ for different gap-junctional coupling strengths $\varepsilon_{g}$ and reiterates that the model cannot match the long wavelength without sufficiently strong gap-junctional coupling. Our neuromechanical model thus relies on gap-junctional coupling between neural modules to achieve the long wavelengths in low-viscosity environments.

A single, isolated neuromechanical module is defined as a neural subcircuit, the corresponding muscles and body section, and local proprioceptive feedback (without coupling through the body or neural circuitry). The dynamics for this isolated module are governed by

Note that this is the model described in Section 2, omitting the intermodular coupling. The isolated modules exhibit robust oscillations over a wide range of parameters, and a single period of the module is shown for each state-variable in Figure 7(a). Thus, the neuromechanical modules are oscillators, wherein each B-class neuron promotes either a dorsal or ventral bend and the local proprioceptive feedback acts to switch the bistable B neurons from one state to the other. The basic cycle of the oscillator is as follows: when activated, the ventral B-class neuron ($V_{V}$) excites the ventral muscles which build up activity ($A_{V}$) to induce a ventral bend (negative $\kappa$); when the curvature $\kappa$ is sufficiently large, the local proprioceptive feedback deactivates the ventral B-class neuron and activates the dorsal B-class neuron, and the cycle continues towards a dorsal bend.

The system of six identical, uncoupled neuromechanical oscillators is described by

where

and $S(\underline{\mathbf{X}})$ is given by equations 25-29. The oscillations correspond to a $T$-periodic limit cycle $\underline{\mathbf{X}}^{LC}(t)$ in $(\kappa,A_{V},A_{D},V_{V},V_{D})$-state-space. This limit cycle can be parametrized by phase

with the initial phase $\theta_{j}^{0}\in[0,1)$. As $\theta_{j}$ increases at a constant rate $\omega=1/T$, $\underline{\mathbf{X}}^{LC}(\theta_{j})$ traces out the limit cycle through state-space and the state of each oscillator on the limit cycle is given by

where the only distinguishing feature between the oscillators is their unique phase $\theta_{j}$. Figure 7(a) shows the components of $\underline{\mathbf{X}}^{LC}(\theta)$.

Rearranging equations 17-23, the neuromechanical model can be written as a network of coupled oscillators:

where $C_{j}(\underline{\mathbf{X}}_{1},\dots,\underline{\mathbf{X}}_{6})$ describes the coupling dynamics from all the modules to the $j^{th}$ module through gap-junctions, nonlocal proprioception, and body mechanics:

The parameter $\varepsilon_{m}=\alpha\mu_{f}\ell^{4}/\mu_{b}$ is the effective mechanical coupling strength.

The oscillations of the isolated module (equations 25-29) are indistinguishable (on the scale of Figure 7) from the oscillations of a module within the fully-coupled network (equations 17-21) at both low and high external fluid viscosity $\mu_{f}$. That is, the dynamics of the coupled module never deviate substantially from the limit cycle of the uncoupled module. This indicates that the intrinsic dynamics of the module dominate over the influence of coupling on any given cycle, which implies that the coupling is “weak”. However, small changes in the phase of a module due to coupling can accumulate over many cycles to significantly influence the phase differences between the modules.

Because the coupling is weak (as defined above), the theory of weakly coupled oscillators can be applied (see [36] for details). The coupling only alters the phase of the oscillators on their respective limit cycles and the effect on amplitude is negligible, therefore the phase completely describes the state of a neuromechanical module. Equation 34 can be reduced to the so-called phase equations, a set of differential equations describing the evolution of the phases of each oscillator:

where $\theta_{j}$ is the phase of the $j^{th}$ oscillator, $\omega$ is the intrinsic frequency, and $H(\phi)$ are the interaction functions that describe the change in frequency (resulting from either mechanical, proprioceptive, or gap-junction coupling) as a function of the phase difference $\phi=\theta_{k}-\theta_{j}$ of a given pair of oscillators:

Here, $Z_{\kappa}(t)$, $Z_{V_{V}}(t)$, $Z_{V_{D}}(t)$ are the $T-$periodic phase response functions to perturbations in the corresponding state variable.

The coupling modalities define the structure of the interaction functions, through the state variables that are coupled, as well as the coupling topology (the connectivity matrices $D_{4}^{-1}$, $W_{g}$, and $W_{p}$ in equation 36). Note that there is a separate H-function for each of the three coupling modalities and these three coupling modalities add linearly to produce the full interaction of the modules. Therefore, the relative contributions of the various coupling types can be analyzed independently through varying the different coupling strengths: fluid viscosity $\mu_{f}$ (through $\varepsilon_{m}$) for mechanical, $\varepsilon_{p}$ for proprioceptive, and $\varepsilon_{g}$ for gap-junctional.

Analyzing a pair of two coupled oscillators gives considerable insight into the coordination that each coupling modality produces separately and the mechanisms of coordination. With only two oscillators, the phase model reduces to

In the two oscillator case, the matrix $D_{4}^{-1}$ is symmetric, so $(D_{4}^{-1})_{12}=(D_{4}^{-1})_{21}=d_{12}$. By defining

and subtracting equation 40 from equation 41, the dynamics of the two oscillator system can be described by a single differential equation for the phase difference between the two oscillators:

where $G_{m}(\phi)=H_{m}(-\phi)-H_{m}(\phi)$, $G_{p}(\phi)=H_{p}(-\phi)$, and $G_{g}(\phi)=H_{g}(-\phi)-H_{g}(\phi)$ are the pair-wise interaction functions, or G-functions of the pair. The stable phase-locked state of the system $\phi^{*}$ is given by $G(\phi^{*})=0,$ $G^{\prime}(\phi^{*})<0$.

We simulate the six-oscillator phase model in order to (i) assess the predictive power of the phase model by a quantitative comparison to the full six-module neuromechanical model and (ii) determine whether the mechanism of gait adaptation analyzed in the two-module case extends to the full six-module case.

Figure 12(a) shows the wavelengths for the six-oscillator phase model (line, circles) and neuromechanical model (crosses) as a function of external fluid viscosity $\mu_{f}$ for $\varepsilon_{p}=0.05$ and $\varepsilon_{g}=0.017$ (these coupling strengths were chosen so that the water-wavelength is approximately 1.5). The wavelengths were computed by equation 54 in Appendix A. The phase model and the neuromechanical model agree quantitatively even at high $\mu_{f}$, where the mechanical coupling strength is several orders of magnitude stronger. Figure 12(b) shows the stable phase differences between neighboring modules in the six-oscillator phase model (lines, circles) and neuromechanical model (crosses) as a function of external fluid viscosity $\mu_{f}$. Again, the phase model and the neuromechanical model are in quantitative agreement. Furthermore, Figure 12(b) shows that increasing fluid viscosity affects the phase-locked states in the six-oscillator case in a similar way as in the two-oscillator case. When neural coupling dominates at low viscosity, the stable phase differences are spread out near 0.9, and as fluid viscosity increases, the mechanical coupling strength increases and the stable phase differences decrease towards antiphase. However the phase differences do not reach pairwise-antiphase as in the two-oscillator case. The large variation between the phase differences across pairs of modules is due to the non-uniformity of coupling matrices $D_{4}^{-1}$, $W_{g},$and $W_{p}$. The modules in the middle receive stronger mechanical coupling than the modules at the boundaries; the boundary modules receive less gap-junctional coupling because they have one fewer neighboring module; and the first module gets zero nonlocal proprioceptive feedback because it has no anterior neighboring module.

The general trend of each phase difference between neighboring modules (decreasing from near-synchrony towards antiphase) underlies the wavelength trend of gait adaptation in Figure 12(a) in both the six-oscillator phase model and the neuromechanical model. Figure 13 shows the curvature kymographs generated by the pairwise phase differences of the phase model (from Figure 12) at three selected fluid viscosities (low $\mu_{f}=1$ mPa s, medium $\mu_{f}=348$ mPa s, and high $\mu_{f}=28$ Pa s). The shortening of the wavelength seen in the curvature kymographs is a direct consequence of the decreasing pairwise phase differences. Thus, the results for the two-oscillator case in Section 4.3 extend to the six-oscillator case: the decrease in wavelength in response to increasing fluid viscosity is the result of the corresponding increase in the relative strength of mechanical coupling, which decreases the phase differences between neighboring modules and yields shorter wavelengths.