A Neuromechanical Model of Multiple Network Rhythmic Pattern Generators for Forward Locomotion in C. elegans

Multiple mechanisms contribute to the generation, propagation, and coordination of the rhythmic patterns necessary for locomotion in Caenorhabditis elegans. Current experiments have focused on two possibilities: pacemaker neurons and stretch-receptor feedback. Here, we focus on whether it is possible that a chain of multiple network rhythmic pattern generators in the ventral nerve cord also contribute to locomotion. We use a simulation model to search for parameters of the anatomically constrained ventral nerve cord circuit that, when embodied and situated, can drive forward locomotion on agar, in the absence of pacemaker neurons or stretch-receptor feedback. Systematic exploration of the space of possible solutions reveals that there are multiple configurations that result in locomotion that is consistent with certain aspects of the kinematics of worm locomotion on agar. Analysis of the best solutions reveals that gap junctions between different classes of motorneurons in the ventral nerve cord can play key roles in coordinating the multiple rhythmic pattern generators.

1. Introduction

Understanding how behavior is generated through the interaction between an organism's brain, its body, and its environment is one of the biggest challenges in neuroscience (Chiel and Beer, 1997; Chiel et al., 2009; Krakauer et al., 2017). Understanding locomotion is particularly critical because it is one of the main ways that organisms use to interact with their environments. Moreover, locomotion represents a quintessential example of how behavior requires the coordination of neural, mechanical, and environmental forces. Caenorhabditis elegans is a particularly ideal candidate organism to study the neuromechanical basis of locomotion because of the small number of neurons in its nervous system and the reconstruction of its neural and muscle anatomy at the cellular level, which has led to a detailed map of the connectivity of its the nervous system (White et al., 1986; Varshney et al., 2011; Cook et al., 2019). However, despite the available anatomical knowledge, how the rhythmic patterns are generated and propagated along the body to produce locomotion is not yet fully understood (Cohen and Sanders, 2014; Gjorgjieva et al., 2014; Zhen and Samuel, 2015).

As with many other organisms, there are likely multiple mechanisms, intrinsic and extrinsic to the nervous system, contributing to the generation, propagation, and coordination of rhythmic patterns that are necessary for locomotion in C. elegans (Rossignol et al., 2006). Until recently, the majority of experimental work on C. elegans locomotion had been focused on understanding the role of extrinsic contributions, specifically the role of stretch-receptor feedback. The proposal that stretch-receptor feedback plays an important role in the generation of movement in the nematode dates back to the reconstruction of the connectome (White et al., 1986). There has since been evidence of mechanically gated channels that modulate C. elegans locomotion (Tavernarakis et al., 1997; Yeon et al., 2018), as well as evidence of a direct relationship between body curvature and neural activity (Wen et al., 2012). However, coordinated rhythmic patterns can also be produced intrinsically, while remaining open to modulation through extrinsic feedback. Intrinsic rhythmic pattern generators are known to be involved in a wide variety of behaviors in a number of different organisms, including insect flight, swimming in molluscs, gut movements in crustaceans, and swimming and respiration in vertebrates (Marder and Bucher, 2001; Goulding, 2009; Arshavsky et al., 2016; Katz, 2016; Minassian et al., 2017; Dasen, 2018). In an intrinsic rhythmic pattern generator, the rhythmic pattern can be generated through the oscillatory properties of pacemaker neurons or it can emerge from the interaction of networks of non-oscillatory neurons (Goulding, 2009). Recent experiments have provided support for the role of intrinsic rhythmic pattern generation in C. elegans locomotion (Fouad et al., 2018; Gao et al., 2018; Xu et al., 2018).

It is increasingly acknowledged that simulation models play an important role in elucidating how brain-body-environment systems produce behavior (Chiel and Beer, 1997; Pfeifer and Bongard, 2006; Ijspeert, 2008; Izquierdo, 2018; Cohen and Denham, 2019). In C. elegans, there has been a surge of theoretical work focused on understanding the neuromechanical basis of locomotion. By taking into consideration the mechanics of the body and its interaction with the environment, several computational models have demonstrated that extrinsic pattern generation alone can result in locomotion (Marom and Shahaf, 2002; Putrenko et al., 2005; Chalasani et al., 2007; Li et al., 2014; Gleeson et al., 2018; Pastore et al., 2018; Fenyves et al., 2020). There have been a handful of models that have considered the potential role of intrinsic rhythmic pattern generators in C. elegans locomotion (Karbowski et al., 2008; Deng and Xu, 2014; Kunert et al., 2014; Denham et al., 2018). Some of these models have considered a circuit capable of intrinsic rhythmic pattern generation in the head motorneurons (Karbowski et al., 2008; Izquierdo and Beer, 2018) or in the command interneurons (Deng and Xu, 2014). Some of the models have imposed a neural activation function in the form of a traveling sine wave to drive a mechanical body to produce movement (Denham et al., 2018; Palyanov et al., 2018). Only a few models have considered the generation of rhythmic patterns from networks of motorneurons in the ventral nerve cord (Kunert et al., 2014; Olivares et al., 2018). In this model we consider both the generation of rhythmic patterns from multiple networks of motorneurons in the ventral nerve cord and the dynamic interaction of these neural patterns with the mechanical body and environment to produce movement.

Given the focus on extrinsic contributions to the generation, propagation, and coordination of rhythmic patterns controlling locomotion in C. elegans, current models have left a major question unanswered: Can multiple network rhythmic pattern generators in the ventral nerve cord coordinate their activity to produce the traveling wave necessary for forward locomotion in the absence of stretch-receptor feedback? And importantly, what are the different possibilities for how this could be accomplished in the worm? In this paper, we coupled multiple repeating neural units in the ventral nerve cord (VNC), whose connectivity was derived from the C. elegans connectome (Haspel and O'Donovan, 2011), to a model of the worm's muscular system (Izquierdo and Beer, 2018) and mechanical body (Boyle et al., 2012), which in turn was situated in a simulated agar environment. In order to examine the feasibility of multiple network rhythmic pattern generators to contribute to forward locomotion within the VNC, the current model deliberately leaves out stretch-receptor feedback, it does not allow for the possibility of motorneurons to be pacemaker neurons, and it does not include neurons outside of the VNC to drive the circuit. We used a real-valued evolutionary algorithm to determine values of the unknown parameters of the neural circuit that optimized the ability of the coupled neuromechanical model worm to match as best as possible the ability to locomote forward on agar. To the degree that this is possible, we will learn something about what components of the worm can recreate movement under these limited conditions. Given the unconstrained nature of many problems in biology (Prinz and Marder, 2004; Gutenkunst et al., 2007), instead of looking for one unique model, we ran multiple evolutionary searches as a way to explore the space of parameter configurations that could lead to the behavior. Each successful search produced a distinct set of parameter values, which led to an ensemble of models that we filtered down to those that were most consistent with the worm's behavior. The properties of this ensemble were then analyzed to identify different possible classes of solutions. Detailed analysis of the operation of the representative exemplars suggests hypotheses for mechanisms that can contribute to the generation and propagation of rhythmic patterns for locomotion in the worm.

2. Methods

In this section, we describe each of the components of the model: the physical environment, the mechanical body, and the neuromuscular system; as well as the optimization technique.

2.1. Model

The neuromechanical model (see Figure 1) integrates the neural unit from Olivares et al. (2018), with the muscular model in Izquierdo and Beer (2018), and the physical model of the body and environment in Boyle et al. (2012).

FIGURE 1

Figure 1. The neuromechanical model integrates a mechanical body, with a muscular system, and a repeating neural unit. (A) Following previous work (Boyle et al., 2012), the mechanical body is modeled in 2D cross-section with a set of variable-width interconnected discrete segments. (B) Each segment consists of two cross-sectional rigid rods (black), damped spring lateral elements (gray), and damped spring diagonal elements (dashed). (C) Following previous work (Izquierdo and Beer, 2018), muscles are modeled as elements that lie along the cuticle that can contract and relax in the dorsoventral plane, staggered along the ventral and dorsal sides of the model worm. Muscle force is distributed across all lateral elements they intersect. Seven repeating neural units innervate three ventral and three dorsal muscles, except the most anterior and the most anterior subunits which innervate four muscles. (D) We modeled the motorneurons in the ventral nerve cord as a network composed of repeating identical subunits. The architecture of each subunit extends previous work (Olivares et al., 2018) and is based on the statistical analysis of the VNC motorneurons given the missing connectome data (Haspel and O'Donovan, 2011). One of seven repeating neural subunits is shown in complete detail. Intraunit connections shown in color; interunit connections shown in black. Chemical synapses shown with solid lines; gap junctions shown with dashed lines.

2.1.1. Environment Model

In the laboratory, C. elegans is typically grown and studied in petri dishes containing a layer of agar gel. The gel is firm, and worms tend to lie on the surface. The experiments in this paper focus on worm locomotion in agar. Given the low Reynolds number physics of C. elegans locomotion, inertial forces can be neglected and the resistive forces of the medium can be well-approximated as a linear drag F = −Cv (Niebur and Erdös, 1991; Boyle, 2009; Cohen and Boyle, 2010; Boyle et al., 2012). The tangential and normal drag coefficients for agar used in this model were taken from those reported in Berri et al. (2009) and used in the model of the body that this work builds on Boyle et al. (2012): $C_{\parallel }=3.2\times 10^{-3}$ kg·s⁻¹ and $C_{\perp }=128\times 10^{-3}$ kg·s⁻¹, respectively (Wallace, 1969; Lighthill, 1976; Niebur and Erdös, 1991; Berri et al., 2009; Boyle, 2009; Boyle et al., 2012).

2.1.2. Body Model

When placed on an agar surface, the worm locomotes by bending only in the dorsal-ventral plane. For this reason, the worm body is modeled in 2D cross-section (Figure 1A). The model of the mechanical body is a reimplementation of the model presented by Boyle et al. (2012). The ~1 mm long continuous body of the worm is divided into discrete segments. The width of the segments change along the length of the body as represented in Figure 1A (for details see Boyle et al., 2012). Each segment is bounded by two cross-sectional rigid rods (Figure 1B). The endpoints of the rods are connected to their neighbors via damped spring lateral elements modeling the stretch resistance of the cuticle. The endpoints of the rods are also connected to the adjacent rods on the opposite side via damped spring diagonal elements modeling the compression resistance of internal pressure. The rest lengths, spring constants and damping constants of the lateral and diagonal elements are taken directly from previous work (Boyle et al., 2012), which in turn estimated them from experiments with anesthetized worms (Sauvage, 2007). The forces from the lateral and diagonal elements are summed at the endpoints of the rods and then the equations of motion are written for the center of mass of each rod. Since each rod has two translational (x, y) and one rotational (ϕ) degrees of freedom, the body model has a total of 3(N_seg+1) degrees of freedom. The current model has N_seg = 50, so a total of 153 degrees of freedom. The full set of expressions for forces, as well as all kinematic and dynamic parameters are identical to those in previous work (Boyle, 2009; Boyle et al., 2012).

2.1.3. Muscle Model

Body wall muscles in the worm are arranged as staggered pairs in four bundles around the body (Waterston, 1988; Altun and Hall, 2009). These muscles can contract and relax in the dorsoventral plane. Following previous work (Izquierdo and Beer, 2018), muscles are modeled as elements that lie along the cuticle (Figure 1C). The force of each muscle is distributed across all lateral elements that they intersect. Because adjacent body wall muscles overlap one another in C. elegans, multiple muscles can exert force on the same lateral elements. Since the model is 2D, we combine right and left bundles into a single set of 24 dorsal and 24 ventral muscles. Muscle forces are modeled as a function of muscle activation and mechanical state using simplified Hill-like force-length and force-velocity properties (Hill, 1938).

Following previous work (Boyle et al., 2012; Izquierdo and Beer, 2018), muscle activation is modeled as a leaky integrator with a characteristic time scale (τ_M = 100ms), which agrees with response times of obliquely striated muscle (Milligan et al., 1997). The muscle activation is represented by the unitless variable $A_m^k$ that evolves according to:

$\begin{array}{ll}\frac{d{A}_m^k}{dt}=\frac{1}{{\tau }_M}\left(I_m^k-A_m^k\right)& (1)\end{array}$

where $A_m^k$ is the total current driving dorsal and ventral (k = {D, V}) muscles along the body (m = 1, …, 24). Also following previous modeling work (Boyle et al., 2012) and experimental evidence that electrical coupling between body wall muscle cells plays only a restricted role for C. elegans body bend propagation (Leifer et al., 2011; Wen et al., 2012), inter-muscle electrical coupling is assumed to be too weak and therefore not included in the model.

2.1.4. Ventral Nerve Cord Circuit

As connectome data is incomplete for the ventral nerve cord (White et al., 1986; Varshney et al., 2011; Cook et al., 2019), we relied on a statistical analysis of the motorneurons in relation to the position of the muscles they innervate to model the repeating neural unit along the VNC (Haspel and O'Donovan, 2011). We modeled the VNC as a neural network composed of seven identical subunits (Figure 1D). The anatomy of the repeating subunit was grounded on previous theoretical work, where we demonstrated that a subset of the components present in the statistically repeating unit found in the dataset were sufficient to generate dorsoventral rhythmic patterns (Olivares et al., 2018). The minimal configuration found in that work included motorneurons: AS, DA, DB, VD, VA, and VB; and chemical synapses: DA→DB, DB→AS, AS→DA, AS→VD, VD→VA, and VD→VB. Given that the subunits need to coordinate their rhythmic patterns with neighboring subunits in order to produce forward locomotion, we added the following connections to adjacent neural subunits found in the statistical analysis of the VNC (Haspel and O'Donovan, 2011): AS⊢⊣VA⁺¹, DA⊢⊣AS⁺¹, VB⊢⊣DB⁺¹, where the superscript +1 indicates that the neuron is part of the posterior subunit. We use this notation to refer to interunit connections only; for intraunit connections we leave the superscript out. The minimal configuration found in previous work (Olivares et al., 2018) did not include motorneuron DD because of the lack of outgoing connections to the rest of the motorneurons within the unit, and therefore its unlikeliness to be involved in the generation of network rhythmic patterns. As the current model involves a neuromuscular system, and DD has neuromuscular junctions that allow it to drive the muscles of the worm, we included it. We also included the connections to and from DD present in the statistical analysis of the VNC (Haspel and O'Donovan, 2011), including intraunit connections: DA→DD, VA→DD, VB→DD, and VD⊢⊣DD; and interunit connections: DB→DD⁺¹, and VA⁺¹→DD.

Dorsal and ventral motorneurons in each unit drive the dorsal and ventral body wall muscles adjacent to them, respectively. The input to the body wall muscles is represented by variable $I_m^k$, such that:

$\begin{array}{ll}{I}_m^k={\displaystyle \sum _{i\in N^k}}{\gamma }_m{q}_i{S}_i& (2)\end{array}$

where k denotes whether the body wall muscle is dorsal or ventral (k = {D, V}) and m denotes the position of the muscle along the body (m = 1, …, 24), the set N^k corresponds to the dorsal/ventral motorneurons, {AS, DA, DB, DD}, {VA, VB, VD} respectively. Following previous work (Boyle et al., 2012), an anterior-posterior gradient in the maximum muscle efficacy is implemented by a linearly (posteriorly) decreasing factor, γ_m = 0.7(1−(((m−1)F)/M)), where γ_m is the efficacy for neuromuscular junctions connecting motorneurons to muscle m, and F is the anterior-posterior gain in muscle contraction. q_i corresponds to the neuromuscular junction strength from motorneuron i. Finally, S_i corresponds to the synaptic output for each motorneuron.

The strengths of the connections in the circuit are unknown. The signs of the connections (i.e., whether they are excitatory or inhibitory) were constrained only for the neuromuscular junctions, but not for the chemical synapses between motorneurons. The AS-, A-, and B-class motorneurons are known to be cholinergic, and therefore excitatory to the muscles they innervate; D-class motorneurons are known to be GABAergic, and therefore inhibitory to the muscles they innervate (McIntire et al., 1993; Rand et al., 1997). We constrained the signs of neuromuscular junctions accordingly. Note, however, we did not constrain the signs of the connections between those motorneurons and other motorneurons in the circuit because these are not known (Marom and Shahaf, 2002; Putrenko et al., 2005; Chalasani et al., 2007; Li et al., 2014; Pastore et al., 2018; Fenyves et al., 2020).

2.1.5. Neural Model

Following studies in C. elegans (Goodman et al., 1998; Mellem et al., 2008) and previous modeling efforts (Izquierdo and Lockery, 2010; Izquierdo and Beer, 2013, 2018), motorneurons were modeled as nodes with simple first order non-linear dynamics (Beer, 1995),

$\begin{array}{ll}{\tau }_i\frac{d{V}_i}{dt}=-V_i+{\displaystyle \sum _{j=1}^N}{w}_{ji}\sigma (V_j+{\theta }_j)+{\displaystyle \sum _{j=1}^N}{g}_{ji}\left(V_j-V_i\right)& (3)\end{array}$

where V_i represents the membrane potential of the i^th neuron relative to its resting potential. The time-constant of the neuron is represented by τ_i. The model assumes chemical synapses release neurotransmitter tonically and that steady-state synaptic activity is a sigmoidal function of presynaptic voltage (Wicks et al., 1996; Lindsay et al., 2011; Kuramochi and Doi, 2017), σ(x) = 1/(1+e^−x). θ_j is a bias term that shifts the range of sensitivity of the output function. The synaptic weight from neuron j to neuron i is represented by w_ji. In line with previous theoretical work (Wicks et al., 1996; Izquierdo and Beer, 2013; Kunert et al., 2017), the electrical synapses were modeled as bidirectional ohmic resistances, with g_ji as the conductance between cell i and j (g_ji > 0). The indices i and j used for the chemical synapses and the gap junctions represent each of the motorneurons in the circuit (AS, DA, DB, VD, VA, and VB) and the specific connectivity between them is given by the neuroanatomy (Figure 1D). Self-connections were included in the chemical synapses term to allow for the functional equivalent of active membrane conductances which have been reported for C. elegans neck muscle motor neurons (Goodman et al., 1998; Mellem et al., 2008). This allows the neural model to reproduce the variety of graded activity that has been described in the free-living nematode C. elegans (Goodman et al., 1998; Mellem et al., 2008; Liu et al., 2009; Lindsay et al., 2011). Specifically, by changing the strength of the self-connection on each neuron, that model neuron can be either smoothly depolarized or hyperpolarized from a tonic resting potential (Mellem et al., 2008), or bistable, with non-linear transitions between a resting potential and a depolarized potential (Goodman et al., 1998).

2.2. Numerical Methods

The model was implemented in C++. The neural model was solved by Forward Euler method of integration with a 0.5ms step. The body model was solved using a Semi-Implicit Backward Euler method with a 0.1ms step.

2.3. Evolutionary Algorithm

All neural circuits described in this article were produced using a simple model of evolution known as a genetic algorithm. The parameters to be searched, such as the weights and signs of the connections, were encoded as a vector of real values. A population of such vectors was maintained. Initially, the vectors in this population were randomly generated. In each generation, the fitness of every individual in the population was evaluated. A new generation of individuals was then produced by applying a set of genetic operators: selection, recombination, and mutation. Once a new population had been constructed, the fitness of each new individual was evaluated, and the entire process repeated.

A naive parameterization of our model would contain around 300 neural parameters. However, it makes little sense to work directly with such a large set of unconstrained parameters. Instead, we assumed that the parameters in each repeating VNC neural unit were identical. Altogether, the model was reduced to a total of 44 free parameters. There are 28 parameters that describe each of the 7 neuron classes and neuromuscular junctions: 7 biases, 7 time-constants, 7 self-connections, and 7 neuromuscular junctions. There are 15 parameters that describe the strength and sign (i.e., excitatory/inhibitory) of the connections: 10 weights for intraunit connections: 9 chemical synapses (AS→DA, AS→VD, DA→DB, DB→AS, VD→VA, VD→VB, DA→DD, VB→DD, VA→DD) and one electrical synapse (VD⊢⊣DD); 5 weights for interunit connections: 2 chemical synapses (VA→DD⁻¹, DB→DD⁺¹) and 3 gap junctions (DA⊢⊣AS⁺¹, VB⊢⊣DB⁺¹, AS⊢⊣VA⁺¹). One additional parameter, F, describes the anterior-posterior gain in muscle contraction.

To evaluate the ability of a configuration of neural parameters to produce locomotion, when embodied and situated, we created a fitness function with two components. The goal of the first component was to make a ventral nerve cord neural unit produce a rhythmic pattern. Specifically, the fitness function required that the B-class motorneurons produce a rhythmic pattern and that the frequency of the rhythmic pattern matched what has been observed for body bending in crawling worms:

$\begin{array}{ll}{F}_1={\displaystyle \prod _{j\in \left\{\text{DB},\text{VB}\right\}}}\left(\frac{2}{A*T}{\displaystyle {\int }_0^T}\left|\frac{d{S}_j}{dt}\right|dt\right)\left(1-\frac{|f_j-f_{\text{a}}|}{f_{\text{a}}}\right)& (4)\end{array}$

where A corresponds to a rhythmic pattern amplitude threshold (A = 0.5), S_j corresponds to the output of the motorneuron, T corresponds to the duration of the simulation, f_j is the frequency of neuron j, and f_a is the frequency of bending in the worm (f_a = 0.44Hz, Cohen et al., 2012). The first part of this equation encourages the circuit to produce a rhythmic pattern by maximizing the rate of change of neural activity in DB and VB. The contribution from this component was capped to a value of 1. The second part of the equation is aimed at matching the frequency of the worm.

The goal of the second component of the fitness function was to make the complete neuromechanical model worm move forward by matching its forward velocity to that of the worm on agar:

$\begin{array}{ll}{F}_2=1-\frac{|\overline{V}-V_{\text{a}}|}{V_{\text{a}}}& (5)\end{array}$

where $\overline{V}$ corresponds to the average velocity of the model worm over the duration of the simulation, and V_a corresponds to the average forward velocity of the worm on agar (V_a = 0.22 mm/s) (Cronin et al., 2005).

3. Results

3.1. Generating an Ensemble of Model Worms That Use Multiple Network Rhythmic Pattern Generators for Forward Locomotion on Agar

As the parameters for the physiological properties of neurons and synapses involved in forward locomotion in C. elegans are largely unknown, we used an evolutionary algorithm to search through the space of parameters for different configurations that could produce forward movement on agar in the absence of stretch-receptors or pacemaker neurons. Because evolutionary runs with the full neuromechanical model and environment are computationally costly, we used an incremental approach. During a first stage, isolated ventral nerve cord neural units were evolved to produce a rhythmic pattern using fitness function F₁ (see section 2) (Figure 2A). Once the isolated neural units could produce rhythmic patterns (fitness > 0.99), they were integrated into the complete neuromechanical model and evolved to move forward in agar using a combined fitness function F₁·F₂ (Figure 2B). We ran 160 evolutionary searches with different random seeds. Of these, 104 (65%) reliably found model configurations that matched the body bending frequency and mean velocity of worms performing forward locomotion on agar (Figure 2C). In other words, the evolutionary search consistently found configurations of the neuroanatomical circuit that could produce forward locomotion through the coordination of multiple network rhythmic pattern generators along the ventral nerve cord.

FIGURE 2

Figure 2. Evolutionary algorithm reliably finds model configurations that match worm forward locomotion. (A) Stage 1: Evolution of the isolated neural unit to match B-class neuron rhythmic patterns and frequency. The best and average fitness of the population are shown in blue and orange, respectively. The average over all evolutionary runs is shown in a solid trajectory and the standard deviation is shown as a lighter shade of the respective color. (B) Stage 2: Evolution of the complete neuromechanical model so as to additionally match the worm's instantaneous velocity. Same color coding as in (A). (C) Distribution of best final fitness across evolutionary runs. Sixty-five percent of the evolutionary searches found a solution with a fitness greater than 0.95 (red bar).

In order to focus on the subset of solutions that resembled forward locomotion in C. elegans most closely, we filtered the set of 104 solutions to those that matched an additional set of locomotion features that were not imposed during the evolutionary search. We applied the following three criteria: (a) Relative role of the different neuron classes in forward locomotion; (b) body curvature; and (c) trajectory curvature. Altogether, 15 solutions fulfilled all three filtering criteria (Figure 3D). We discuss each of the criteria in turn.

FIGURE 3

Figure 3. Filtering solutions to include only those that match the ranges reported in the literature about worm locomotion on agar. (A) Relative role of the different neuron classes in forward locomotion. Proportion of the forward locomotion speed maintained by each model worm when the neuromuscular junction from the B-class motorneurons (x-axis) or the A-class motorneurons (y-axis) are ablated. See F₂ in section 2 for the measure of locomotion performance. Each point in the figure represents a single solution. Solutions in the shaded region represent those that match the filtering criteria: (1) Ablation to A-class neuromuscular junctions should not impair forward locomotion entirely; and (2) Ablation to B-class neuromuscular junctions should impair forward locomotion performance. (B) Body curvature. Measures of the model worms' body wavelength (x-axis) and their anterior-posterior curvature profile (y-axis) (see Supplementary Material 1 for details). Green shaded areas represent biologically plausible ranges. In the previous two figures, solutions in the darkest shaded region represent those that match both criteria. Histograms are shown for the criteria on each axis. (C) Trajectory curvature. Distribution of the radius of curvature for the trajectories of each model worm in the 2D plane (see Supplementary Material 1 for details). The blue shaded area represents solutions with a relatively straight trajectory. (D) Venn diagram representing the distribution of the 104 selected solutions according to the fulfillment of the three different filters: relative role of different neural classes in red, body curvature in green, and trajectory curvature in blue. We focus our analysis on the 15 solutions that matched all criteria.

3.1.1. Relative Role of the Different Neuron Classes in Forward Locomotion

A- and B-class neurons have been implicated in backward and forward locomotion, respectively, through ablations performed at the larval stage, when only DA and DB neurons are present (Chalfie et al., 1985). Specifically, these studies have revealed that ablating B-class motorneurons prevents forward locomotion but not backward, and that ablating A-class motorneurons prevents backward but not forward locomotion (Chalfie et al., 1985). More recently, neural imaging studies in the adult have provided evidence that both A- and B-class motorneurons are active during locomotion (Faumont et al., 2011). There is also evidence that the activity of B-class motorneurons is higher during forward locomotion than the activity of A-class motorneurons, and vice-versa during backward locomotion (Kawano et al., 2011). In all evolved solutions of our model, both A- and B-class motorneurons are actively involved in forward locomotion. This is the case because the solutions are all network oscillators. So although we cannot ablate A- or B-class neurons without disrupting the network oscillation, we can silence each of their contributions to the muscles. In order to focus only on the solutions where the B-class input to muscles is necessary to produce forward locomotion but not the A-class, we simulated each solution while eliminating the neuromuscular junctions from B-class motorneurons and from A-class motorneurons, independently. We then evaluated the velocities of the model worms as a result of this manipulation (Figure 3A). We selected solutions that met the following two criteria: (1) eliminating the A-class neuromuscular junction does not seriously compromise locomotion (i.e., velocity greater than 20% of target velocity); and (2) eliminating the B-class neuromuscular junction does compromise forward locomotion (i.e., velocity less than 20% of target velocity). A total of 74 solutions fulfilled both criteria.

3.1.2. Body Curvature

In addition to the frequency of the body bends, there are a number of other features of the kinematics of movement during forward locomotion that have been characterized (Cronin et al., 2005; Karbowski et al., 2008; Pierce-Shimomura et al., 2008; Berri et al., 2009; Fang-Yen et al., 2010; Vidal-Gadea et al., 2011; Cohen et al., 2012; Lebois et al., 2012; Wen et al., 2012; Yemini et al., 2013; Butler et al., 2015; Xu et al., 2018). We further filtered our solutions based on two features: the body-bending wavelength, and the anterior-posterior curvature profile. Measurements of the wavelength of the body during locomotion in agar fall in the range of 0.4 to 0.9 body length (Cronin et al., 2005; Karbowski et al., 2008; Pierce-Shimomura et al., 2008; Berri et al., 2009; Fang-Yen et al., 2010; Vidal-Gadea et al., 2011; Cohen et al., 2012; Lebois et al., 2012; Yemini et al., 2013; Butler et al., 2015). We evaluated the body wavelength in all solutions and selected those that fell within the observe range (Figure 3B). The anterior-posterior curvature profile corresponds to the relative amount of curvature along the body axis and has been shown to be more pronounced near the head of the worm than the tail (Boyle et al., 2012; Wen et al., 2012; Xu et al., 2018). We evaluated the mean curvature in the anterior-posterior axis in all solutions and selected those with a negative slope in the linear regression that fit the curvature profile (Figure 3B). Altogether, we narrowed down the 104 solutions to 30 that fulfilled both criteria (Figure 3D).

3.1.3. Trajectory Curvature

The translational direction of C. elegans during forward locomotion tends to be relatively straight, with only a small degree of curvature in the absence of stimuli (McIntire et al., 1993; Peliti et al., 2013). In the evolved model worms, the straightness in the trajectory was not optimized, so the distribution of curvature in the translational trajectory is broad (Figure 3C). In order to filter out model worms that curved much more than the worm during forward locomotion, we measured the radius of curvature for the trajectories of the centers of mass of each model worm in the 2D plane (see Supplementary Material 1 for details). We set a threshold of 1 mm in trajectory curvature radius (Figure 3C) and we found 77 solutions that moved as straight as the worm (Figure 3D), even in the absence of proprioceptive information.

3.2. Behavior of Model Worms

Despite the absence of stretch-receptor feedback and pacemaker neurons, when simulated, all 15 selected model worms exhibited regular dorsoventral bends that propagated from head to tail (see example of one model worm in Figure 4). The kymograph shows that the traveling wave is not perfectly smooth across the body, instead there is some amount of punctuation. This is due to a combination of simplifying factors, including that the model incorporates seven distinct neural subunits and the lack of proprioceptive feedback. It can also be seen that the traveling wave is less pronounced in the head and tail regions of the body. This is because the mechanical body itself smooths the distinct rhythmic patterns by each neural unit in relation to the unit anterior and posterior to it, which is not the case for the first and last units. The locomotion behavior of all 15 solutions can be observed in animations provided in the Supplementary Material 2. A number of recent experiments have provided support for the possibility that multiple intrinsic rhythmic pattern generators in the ventral nerve cord could be involved in aspects of forward locomotion in the worm (Fouad et al., 2018; Xu et al., 2018). In this section, we examine the different ways in which these model worms are consistent with some aspects of what has been observed in those experiments.

FIGURE 4

Figure 4. Locomotive behavior of model worms. Example from one model worm from the ensemble of solutions. (A) Kymograph depicting dorsoventral bending across the body (y-axis) and over time (x-axis). The intensity of red and blue depict dorsal and ventral curvature, respectively (see Supplementary Material 1 for the method used to calculate curvature). (B) Instantaneous velocity of the body center of mass (black trace) in relation to average velocity used for the fitness function (red line). (C) Neural activity for dorsal neurons in the first neuromuscular segment. AS in black, DA in red, DB in blue and DD in green. (D) Neural activity for ventral neurons in the first neuromuscular segment, VA in red, VB in blue and VD in green. (E) Worm body posture at different points in time (in seconds).

3.2.1. Posterior Rhythmic Patterns Persist Despite Anterior Paralysis

Recent experiments have provided evidence that posterior dorsoventral bending persists despite anterior paralysis (Figure 2 in Fouad et al., 2018; and Figure 3 and Supplementary Figure 3 in Xu et al., 2018). The model presented here is consistent with their experimental finding. In order to demonstrate this, we replicated the experimental condition on the model worms in two ways. First, we suppressed the neuromuscular junction activity for the three anterior-most neural units. Second, we silenced the neural activity of all neurons in those same three anterior-most neural units. Note that taking into consideration the suppression of stretch-receptor feedback was not necessary given that this model did not include stretch-receptor feedback. We examined the resulting kinematics of movement under both conditions. Specifically, we measured the magnitude of the amplitude of dorsoventral rhythmic patterns in the head and the tail. In both conditions, we observed a sharp reduction in dorsoventral bending in the head, but only a slight reduction of dorsoventral bending in the posterior regions of the body in all 15 solutions (Figure 5A). Furthermore, coordination of the multiple rhythmic patterns in the posterior part of the body remained intact (see example from one model worm in Figure 5B). Therefore, as with the worm, posterior dorsoventral bending persists in model worms despite anterior paralysis.

FIGURE 5

Figure 5. Filtered models are consistent with recent experimental observations. (A) Rhythmic posterior undulation persists despite anterior paralysis. Total bending amplitude (y-axis) evaluated as in previous work (Fouad et al., 2018). [i] When motorneurons in the head are ablated, bending amplitude in the head decreases but not in the tail. [ii] Similarly, when body-wall muscles (BWM) in the head are inactivated, bending amplitude in the head decreases but not in the tail. Therefore, paralysis in the head does not abolish bending in the tail. (B) Example kymograph from one model worm shows bending over time when neuromuscular junctions in the head are inactivated. The intensity of red and blue depict dorsal and ventral curvature, respectively. The color-coding is the same as the one used in Figure 4A. (C) Rhythmic undulation persists simultaneously in the head and the tail despite midbody paralysis. [i] When motorneurons in midbody are ablated, bending amplitude decreases in the midbody but not in the head or the tail. [ii] Similarly, when body-wall muscles (BWM) in the midbody are inactivated, bending amplitude decreases in the midbody but not in the head or the tail. Therefore, midbody paralysis demonstrates that head and tail are capable of simultaneous and uncoordinated rhythmic patterns. (D) Example kymograph from one model worm shows bending over time when neuromuscular junctions in the midbody are inactivated. The color-coding is the same as the one used in Figure 4A. (E) Overexpression of electrical synapses on B-class motorneurons induce complete body paralysis. Overexpression was simulated by increasing the synaptic strength of B-class gap junctions: VB⊢⊣DB⁺¹, DB⊢⊣DB⁺¹, VB⊢⊣VB⁺¹. (F) Speed as a function of gap junction overexpression in the model worms. The disks represent the mean and the bars represent the standard deviation. As the synaptic strength of the B-class gap junctions is increased, the speed of the model worms decreases.

3.2.2. Head and Tail Are Capable of Simultaneous and Uncoupled Rhythmic Patterns

Recent experiments have also provided evidence that the head and the tail are capable of simultaneously producing uncoupled rhythmic patterns (Figure 3 in Fouad et al., 2018). The model presented here is also consistent with this key component of their experimental finding. To demonstrate this, we again replicated the experimental condition in two ways. First, we suppressed the neuromuscular junction activity for the three mid-body neural units. Second, we silenced the neural activity of all neurons in those same three mid-body neural units. In both conditions, we observed that suppressing mid-body components did not eliminate body bending in either the head or the tail (Figure 5C). In other words, like in the worm, uncoupled dorsoventral rhythmic patterns were present simultaneously in the head and the tail (see example from one model worm in Figure 5D). It is important to note that in the experiments on the worm (Fouad et al., 2018; Xu et al., 2018), the frequency of oscillations in the anterior and posterior were different. By design, this aspect of their findings cannot be replicated by our model because every unit is identical and because neurons outside the VNC were not included. In additional experiments, we enabled and disabled each of the segments individually. We found that a single segment is not sufficient to drive bending and locomotion. Also, other than the most posterior segments, most segments contribute to forward locomotion (see Supplementary Material 6).

3.2.3. Strengthening Gap Junctions in B-Class Motorneurons Impairs Locomotion

Finally, recent experiments have provided evidence that strengthening the gap junctions (via genetic overexpression of unc-9, one of the genes responsible for gap junctions in C. elegans) in the B-class motorneurons leads to constitutive paralysis in the worm (Xu et al., 2018). Although their experiments involved the manipulation of both local electrical coupling between motorneurons as well as the descending electrical coupling from the command interneurons, the authors of that work suggest that the strong electrical couplings between the motorneurons tends to synchronize motor activity along the whole body, thus deteriorating the ability to propagate the bending wave. Our model is consistent with their interpretation of their results. To demonstrate this, we systematically increased the synaptic weight of the gap junctions interconnecting B-class motorneurons (both the VB⊢⊣DB⁺¹, VB⊢⊣VB⁺¹, and DB⊢⊣DB⁺¹) and measured the resulting bending along the body and speed of the model worm. As the strength of the gap junctions was increased, the bending in the body decreased. Noticeably, the effect was more pronounced in the tail than in the head (Figure 5E). Accordingly, the velocity of the simulated worms also decreased as the strength of the gap junctions were increased, leading ultimately to total lack of movement forward (Figure 5F). The reason for the reduced velocity is the increased synchronization of the different network oscillators along the body caused by the increased strength of the gap junctions between them. Therefore, in the model worms, strengthening the electrical couplings between the motorneurons deteriorated their ability to propagate the bending wave. There are some parallels between this result and the experiments described previously (Xu et al., 2018). However, it is important to keep in mind the limitations of this match. First, in those experiments direct testing of the functional contribution of local electrical couplings was not possible because of experimental limitations, and therefore the electrical coupling between AVB and the B-class motorneurons could also be playing a role. Second, in those experiments stretch-receptor feedback has not been eliminated, and it could therefore also be playing a role.

3.3. Rhythmic Pattern Generation and Coordination in the Ensemble of Model Worms

In the previous section, we provided evidence that the simulated model worms can produce locomotion without stretch receptors and without pacemaker neurons in a way that both resembles the kinematic characterization of the worm's forward movement on agar and is consistent with various experimental manipulations. This suggests that the way these model worms operate could be illustrative for understanding the mechanisms responsible for locomotion in the worm. Two basic mechanisms are necessary for a chain of multiple rhythmic pattern generators to drive locomotion in the worm. First, a network of neurons must be able to generate rhythmic patterns intrinsically. Second, adjacent rhythmic pattern generators must coordinate their activity with the appropriate phase delay along the anterior-posterior axis. In what follows, we examine the model worms in detail to answer the following three questions: How dominant are inhibitory or excitatory connections in the evolved rhythmic pattern generators? How do the model worms generate dorsoventral rhythmic patterns? And how do they coordinate these rhythmic patterns across the length of the body to generate a propagating wave capable of producing thrust against an agar surface?

3.3.1. Excitatory-Inhibitory Pattern in the Ensemble of Model Worms

It has been proposed that generating intrinsic network rhythmic patterns is difficult because the network would have to rely extensively on inhibitory connections (Cohen and Denham, 2019). However, the evolutionary search revealed multiple instantiations of possibilities over a wide range of the inhibition/excitation spectrum (see Supplementary Material 3 for the full set of parameters for each of the solutions). Seven out of the 15 solutions contained a majority of excitatory synapses. Furthermore, across the 15 models analyzed, it is possible to find one with any of its six chemical synapses in an excitatory configuration. This suggests a wide range of possibilities for the feasibility of multiple intrinsic network rhythmic pattern generators in the ventral nerve cord.

3.3.2. Model Worms Use the Dorsal AS-DA-DB Subcircuit to Generate Rhythmic Patterns

How do these model worms generate rhythmic patterns? To answer this question, we first determined which set of neurons are involved in producing rhythmic patterns. For a subcircuit to be capable of generating rhythmic patterns in the absence of pacemaker neurons, a recurrently connected set of neurons are required. There are three possible subcircuits in the VNC unit that are capable of generating intrinsic network rhythmic patterns: AS-DA-DB, VD-VA-DD, and VD-VB-DD (Figure 6A). We examined whether each of these subcircuits alone could produce rhythmic patterns (Figure 6B). We measured the total cumulative change in neural activity as an indicator of rhythmic patterns. Because the subcircuits were evaluated in isolation from the rest of the network, we examined each of them with a wide range of compensatory tonic input to each neuron. Consistent with our previous work in the isolated neural unit (Olivares et al., 2018), all model worms generated rhythmic patterns in the AS-DA-DB subcircuit. Out of the 15 model worms examined, only one of them (solution M11) generated rhythmic patterns also in the VD-VB-DD subcircuit and only one (solution M6) generated weak rhythmic patterns in the VD-VA-DD subcircuit.

FIGURE 6

Figure 6. Rhythmic patterns originate primarily in the dorsal core subcircuit: AS-DA-DB. (A) Subcircuits within a neural unit where network pattern generation is possible. (B) Ability to produce rhythmic patterns for the three subcircuits in each of the solutions from the ensemble. The ability to produce a rhythmic pattern was estimated using the average of the absolute value of the derivative of the neural outputs over time, normalized to run between 0 and 1.

3.3.3. Model Worms Use the AS→VD Connection to Propagate the Rhythmic Patterns to the Ventral Motorneurons

Despite the primary role of the dorsal motorneurons in the generation of rhythmic patterns, all model worms show rhythmic patterns in their ventral neural traces. In the majority of model worms (11 out of 15), the ventral motorneurons in an isolated subunit can produce rhythmic patterns (Figure 7A). How does the rhythmic pattern propagate to the ventral motorneurons in these model worms? There are two possibilities: the AS→VD or the DA→DD chemical synapses. We examined whether the ventral B-class motorneuron could produce rhythmic patterns when either of those connections was ablated (Figure 7). In 10 out of the 11 solutions the AS→VD (and not the DA→DD) synapse was necessary to propagate the rhythmic patterns from the dorsal core (Figures 7B,C). This is consistent with our previous work in the isolated neural unit (Olivares et al., 2018). In one of the model worms (solution M11), neither of the connections were necessary. Recall from the previous section that this solution was the only one that also generated strong rhythmic patterns in the ventral side. There are four solutions where the ventral motorneurons do not exhibit a rhythmic pattern in the isolated subunit. In these solutions, the rhythmic pattern in the ventral motorneurons is due to interunit contributions.

FIGURE 7

Figure 7. Model worms use the AS→VD connection to propagate the rhythmic pattern to the ventral motorneurons. (A) Rhythmic patterns in VB motorneurons. Four solutions do not exhibit rhythmic patterns in VB (solid black). (B) Rhythmic patterns in VB is abolished when the connection AS→VD is ablated (connection shown in red). This is true in all solutions, except M11 (shown with an asterisk), which shows rhythmic patterns in the ventral subcircuit. (C) Rhythmic patterns in VB persists when the connection DA→DD is ablated (connection shown in red).

3.3.4. Rhythmic Patterns Coordinate Through a Combination of Three Key Interunit Gap Junctions: AS-VA⁺¹, DA-AS⁺¹, and VB-DB⁺¹

That the model worms move forward is evidence that the multiple rhythmic pattern generators along the body coordinate their activity. But how is the coordination between the different units achieved? To answer this question, we examined the necessity and sufficiency of each interunit connection, chemical and electrical (Figure 8). In all the model worms examined, only the gap junctions played a role in coordinating rhythmic patterns among the different units in the VNC. The interunit chemical synapses were neither necessary nor sufficient for the coordination. In 9 of the 15 solutions examined, a single gap junction was both necessary and sufficient to coordinate the chain of multiple rhythmic pattern generators to drive locomotion forward in the worm. The VD⊢⊣DB⁺¹ gap junction was necessary and sufficient to coordinate rhythmic patterns in four of the solutions; the DA⊢⊣AS⁺¹ gap junction was necessary and sufficient in three solutions; the AS⊢⊣VA⁺¹ gap junction was necessary and sufficient in two. These solutions are particularly interesting because of how simple they are (Figure 8). We analyze one from each of these groups in more detail in the next section. There are three solutions where multiple single gap junctions are sufficient, but no single gap junction was necessary. These solutions use redundant mechanisms to coordinate (Figure 8). Finally, there are three solutions where no single connection is sufficient but several of them are necessary. These solutions are the most complex of the ensemble because they rely on multiple gap junctions to coordinate (Figure 8).

FIGURE 8

Figure 8. Necessity and sufficiency of each interunit connection to coordinate rhythmic patterns for locomotion in each of the model solutions. Each of the 15 solutions (labeled M1-M15) are shown on the x-axis. On the y-axis is the locomotion performance of each solution as a result of examining each of their interunit connections (labeled by color, see legend) for necessity and sufficiency. Locomotion performance was measured as the ability of model worms to match the speed of the worm (see F₁ in section 2). To evaluate necessity, we ablated the connection in question and examined the worm's ability to move forward (solid disks). To evaluate sufficiency, we ablated all but the connection in question and again examined the worm's ability to move forward (circles). Analysis of all 15 solutions revealed three categories of strategies for coordination. “Simple” solutions correspond to those in which a single gap junction is both necessary and sufficient to coordinate the chain of multiple rhythmic pattern generators that drive locomotion. These group of solutions are further subdivided based on which of the three gap junctions is responsible for coordinating the subunits: VD⊢⊣DB⁺¹, DA⊢⊣AS⁺¹, and AS⊢⊣VA⁺¹. “Redundant” solutions are those in which more than one solution is sufficient to drive locomotion. “Complex” solutions are those in which no single gap junction is responsible for coordinating between units. Asterisks in the x-axis label mark the solutions with the highest single sufficient connection from each of the solutions in the “Simple” groups.

3.4. Analysis of Individual Representative Solutions

We have analyzed the properties of the ensemble and we have identified different possible categories of solutions based on how they coordinate the multiple rhythmic patterns. In order to understand how the circuits in the ensemble work, we need to move away from the general features of the ensemble and instead analyze in detail the operation of specific circuits. We selected three representative solutions from the “simple” group to analyze in detail, one belonging to each different cluster of solutions based on which gap junction was responsible for coordinating the subunits. Individuals were selected based on the highest performance on the sufficiency test (i.e., solution M14 for VD⊢⊣DB⁺¹, solution M15 for DA⊢⊣AS⁺¹, and solution M6 for AS⊢⊣VA⁺¹). Based on the results from the previous section, we simplified solutions to their minimal circuit configurations. Simulated models could still perform locomotion efficiently in these simplified configurations (Figure 9A). In all three simplified solutions, the kinematics of movement exhibit a rhythmic pattern in the head that travels posteriorly in a way that remains consistent with what has been observed in the worm (Figure 9B). Because all three solutions can generate movement forward, we know that the multiple rhythmic pattern generators along the body coordinate to achieve the required phase shift. From the previous section we also know that an individual synapse is sufficient to coordinate the rhythmic patterns. In this section, we examine how the coordinated phase-shift is achieved in each of these solutions.

FIGURE 9

Figure 9. Mechanisms of anterior-posterior coordination. (A) Minimal network capable of driving locomotion in each of the solutions from the “simple” group: M14 for VD⊢⊣DB⁺¹ (A1), M15 for DA⊢⊣AS⁺¹ (A2), and M6 for AS⊢⊣VA⁺¹) (A3). Arrows represent excitatory chemical synapses. Connections ending in circles represent inhibitory chemical synapses. Connections with line endings represent gap junctions. (B) Kymographs for each of the minimal configurations above show coordinated bending waves through the body. The intensity of red and blue depict dorsal and ventral curvature, respectively. The color-coding is the same as the one used in Figure 4A. (C) Entrainment analysis for each of the solutions reveals the directionality of the coordination among the subunit rhythmic pattern generators. The purple trajectory depicts the shift in phase that occurs in the posterior-most unit when the phase of the anterior-most unit is displaced. The brown trajectory depicts the shift in phase that occurs in the anterior-most unit when the phase of the posterior-most unit is displaced. In solutions M14 and M15, the anterior-most neural unit is capable of entraining the posterior-most neural unit but not the other way around (C1,C2). This suggests the coordination afforded by these two gap junctions is directed posteriorly. On the contrary, in solution M6, it is the posterior-most neural unit that can entrain the anterior-most unit and not the other way around (C3). This suggests the coordination afforded by this gap junction is directed anteriorly.

3.4.1. Directionality of Coordination

The first thing we need to understand about coordination in these circuits is their directionality. Do anterior units influence the ones posterior to them, or vice-versa? Because the neural units along the VNC are coordinating their phases through gap junctions that allow for bi-directional communication, the directionality of coordination is not directly obvious.

First, in the solution that relies on the VB⊢⊣DB⁺¹ gap junction (Figure 9A1), the anatomy suggests that the rhythmic pattern propagates posteriorly. This is because the interunit connection VD⊢⊣DB⁺¹ places the posterior neural subunit effectively downstream of the anterior neural subunit. The rhythmic pattern in the anterior dorsal core propagates ventrally. Then the VB⊢⊣DB⁺¹ gap junction coordinates the rhythmic pattern with the dorsal core unit immediately posterior to it. Therefore, in this solution, despite the bi-directionality of the coordinating gap junction, the anterior units are likely to be setting the phase of the posterior ones, and not the other way around. In order to test this hypothesis, we performed an entrainment analysis. We introduced a shift in phase first in the anterior-most neural unit and then in the posterior-most neural unit, and we measured the degree to which the rest of the neural units adopted the new phase (Figure 9C1). As expected, when the phase was shifted in the anterior-most unit, the rest of the body adopted that shift successfully; when the phase was shifted in the posterior-most unit, the rest of the body was unaffected.

Second, in the solution that relies on the AS⊢⊣VA⁺¹ gap junction (Figure 9A3), the anatomy suggests that the rhythmic pattern propagates anteriorly. This change in directionality is a result of the interunit connection AS⊢⊣VA⁺¹, placing the anterior rhythmic pattern generator downstream of the posterior one. The rhythmic pattern in the posterior dorsal core, once propagated ventrally, affects through the AS⊢⊣VA⁺¹ gap junction the rhythmic pattern of the dorsal core in the unit immediately anterior to it. Therefore, opposite to the previous model worm, in this model worm the posterior units are likely to be setting the phase of the anterior ones, and not the other way around. This is again despite the bi-directionality of the coordinating gap junction. We tested this hypothesis using the same entrainment analysis as before (Figure 9C3). As expected, when the phase was shifted in the anterior-most unit, the rest of the body was unaffected; when the phase was shifted in the posterior-most unit, the rest of the body adopted that shift successfully.

Finally, in the solution that relies on the DA⊢⊣AS⁺¹ gap junction (Figure 9A2), anatomy alone cannot tell us whether the coordination is occurring anteriorly or posteriorly. Because the connection is directly between the two neural subunits (neither one is downstream of the other), and because the coordinating component is a bi-directional gap junction, the coordination can occur in either direction. We used the same entrainment analysis as before to examine the directionality of coordination in this model worm (Figure 9C2). When the phase was shifted in the anterior-most unit, the rest of the body adopted that shift successfully; when the phase was shifted in the posterior-most unit, the rest of the body was unaffected. Thus, in this model worm the coordination of the shift occurs from head to tail.

3.4.2. Interunit Phase-Shift

The second aspect of the coordination that is crucial to understanding how locomotion is generated is the shift in phase between adjacent neural units. In order to examine this, we identified the approximate shift in phase that occurs at every step of the way from the DB neuron in one unit to the DB neuron in the adjacent unit (Figure 10). We selected to measure the shift in phase between adjacent B-class neurons because of their primary role in forward locomotion. Although the neural dynamics in the model correspond to periodic patterns of activity, the specific shape of each neural activity is different. Because of the differences, we cannot relate the dynamic of two neurons as merely a shift in phase [i.e., f(t) = g(t+T), where f and g are the dynamics of the two neurons]. Nevertheless, we can approximate the shift in phase by assuming that the neurons in the model share the same rhythmic pattern frequency. This is the case particularly in the midbody subunits. For this analysis, we used units 3 and 4 to calculate the shift in phase. In order to estimate the phase of each neuron, we calculate the middle point between the maximum and minimum rate of change for one rhythmic cycle for each neuron.

FIGURE 10

Figure 10. Phase delay among adjacent subunits in the three representative solutions analyzed in more detail: M14 for VD⊢⊣DB⁺¹ (A), M15 for DA⊢⊣AS⁺¹ (B), and M6 for AS⊢⊣VA⁺¹ (C). The highlighted neurons and connections for each network illustrates the shortest path from one DB to the DB in the adjacent neural unit in the direction of the transmission of information determined from the directionality analysis in the previous section. Despite large differences in how these different model circuits operate, the shift in phase from one unit to the one immediately posterior to it are relatively similar.

The highlighted neurons and connections for each network illustrates the shortest path from one DB to the DB in the adjacent neural unit in the direction of the transmission of information determined from the directionality analysis in the previous section. The first thing to note is that the path is different for each network (Figure 10). Second, the shift in phase between two neurons is different among the different networks. For example, in one of the solutions (Figure 10A), DB→AS is linked to an eight-degree shift in phase between DB and AS, whereas the same connection is linked to a 72-degree shift in phase in another network (Figure 10B). However, despite the differences in how the model circuits operate at the level of pairwise neuron interactions, the shifts in phase from one complete neural unit to the one immediately posterior are relatively similar. It is key to note that it is this shift in phase from one unit to the next which is the primary functional activity of the network, which ultimately leads to efficient forward locomotion on agar. This analysis suggests that high level of variability at the neuron-level implementation of solutions can result in similar functional results. To highlight this result, we examined the rest of the circuits in the “simple” networks. The phase shift measured between neurons DB→AS had a mean of 154.0 degrees and a standard deviation of 109.9. Yet, adjacent units had a mean phase shift of 54.4 with a standard deviation of only 8.9.

3.5. Interunit Gap Junctions Present in Connectome

It is important to recall that the ventral nerve cord connectome data is to date still incomplete (White et al., 1986; Varshney et al., 2011; Cook et al., 2019). Therefore, the repeating neural unit upon which we based our model is based on a statistical summary of the VNC (Haspel and O'Donovan, 2011). In this section, we address how the key components that we have identified map onto the actual neuroanatomy of the worm. We examined the most recent reconstructions of the hermaphrodite and the male (Varshney et al., 2011; Jarrell et al., 2012; Xu et al., 2013; Cook et al., 2019) for the existence of the three key interunit gap junctions responsible for coordinating the multiple rhythmic patterns in the model worms. We found that all three key components occur in a large portion of the ventral nerve cord in both hermaphrodites and males (Figure 11). Moreover, because the connectome reconstruction is still incomplete in the mid-body and posterior section of the VNC, additional connections are likely to be present in these regions. The presence of additional connections can only strengthen the possibility of network oscillations. Although no single connection is present for every pair of adjacent units, there are combinations of different connections scattered throughout the full length of the ventral nerve cord. As observed in many of the evolved solutions, these three interunit gap junctions can work together to help coordinate rhythmic patterns. One thing to keep in mind is that not all of these connections are directed in the way we have idealized in the model. For example, although most AS⊢⊣VA connections are directed posteriorly (i.e., with AS anterior to VA), in one of the connections the directionality is inverted (see posterior of the hermaphrodite). Assuming the directionality shown in the model, VB⊢⊣DB and DA⊢⊣AS will coordinate units posteriorly, and AS⊢⊣VA will coordinate units anteriorly. Variability in the directionality in the connectome will have the effect of inverting the directionality of the coordination for those connections. Altogether, from these results, it seems plausible that multiple network rhythmic pattern generators in the ventral nerve cord could robustly coordinate their phases, even in the absence of stretch-receptor information, both anteriorly and posteriorly in the worm.

FIGURE 11

Figure 11. Key interunit gap junctions present in the C. elegans connectome for the hermaphrodite (A) and the male (B). Gap junctions are color coded based on the pair of motorneurons they connect, as in previous figures: VB⊢⊣DB in yellow, DA⊢⊣AS in green, and AS⊢⊣VA in red. Connections shown are present in one or several of the main datasets. For the hermaphrodite, each connection is marked as appearing in one or more of the following three datasets: j (Jarrell et al., 2012), v (Varshney et al., 2011), or c (Cook et al., 2019). For the male, each connection is marked as appearing in one or more of the following two datasets: j (Jarrell et al., 2012), or c (Cook et al., 2019).

4. Discussion

Given the importance of locomotion for the worm, there is a likelihood that the worm has multiple, redundant, and overlapping mechanisms for generating and coordinating rhythmic patterns. One of the mechanisms for which there is a growing consensus is proprioception (Tavernarakis et al., 1997; Mailler et al., 2010; Boyle et al., 2012; Cohen et al., 2012; Wen et al., 2012; Fieseler et al., 2018; Izquierdo and Beer, 2018). In this paper, we examined the possibility that multiple intrinsic network rhythmic pattern generators could coordinate through motorneuron gap junctions to produce a traveling wave along the body. The model did not aim to reproduce any one existing experiment in the literature. Instead, the model aimed at proposing a new experiment that has not been possible experimentally yet: to eliminate proprioception in the worm, to eliminate descending inputs from command interneurons, to eliminate the ability for any motorneurons to produce pacemaking activity, and to examine the worm's ability to propagate a bending wave. The hypothesis that this model makes is that it is theoretically possible for the generation and coordination of oscillatory patterns that can generate a traveling wave, even under those conditions. If this mechanism were to be operational in the worm, then it is likely that it would be an additional contributing mechanism for forward locomotion in the worm.

We developed a model of a repeating ventral nerve cord neural circuit, and embodied and situated it to produce locomotion. In order to test whether coupling between motorneurons could be capable of generating and coordinating rhythmic patterns, the model deliberately excluded proprioceptive feedback, it excluded the ability of A- and B-type motorneurons to exhibit intrinsic oscillatory activity, and it excluded descending signals from the command interneurons. Using an evolutionary algorithm, we found 15 model configurations that reproduced the kinematics of locomotion on agar. The models were consistent with some of the recent experimental results that suggest the feasibility of multiple sources of rhythmic pattern generation (Fouad et al., 2018; Xu et al., 2018).

The multiple intrinsic oscillatory mechanism found by evolution in this model depends entirely on the collective operation of the following ventral nerve cord motor components: It depends on motorneurons AS, DA, and DB to produce the network oscillations; it depends on motorneurons VD, VA, and VB to drive ventral oscillations out of phase to the dorsal ones; and it depends on one of a set of three different inter-segmental gap junctions VB-DB, DA-AS, and AS-VA to coordinate the multiple oscillations across the length of the body. The model provides theoretical support that these components would be key places to look for intrinsic pattern generation and coordination in the worm, if it were to be identified experimentally. It follows naturally that the mechanism proposed is not robust to ablations to any of those individual components. To the degree that ablations to any of those individuals components disrupts the generation and coordination of rhythmic patterns in the worm, this mechanism for generating and coordinating network oscillations is not the sole or primary mechanism responsible for forward locomotion in the worm. In the absence of any one of the VNC motorneurons neurons, this mechanism would cease to be a possible contributor to locomotion. The hypothesis that this model makes is that in the absence of stretch receptor feedback, but otherwise with all neurons intact, it is theoretically possible for the worm to generate and coordinate multiple intrinsic oscillations. This has not been tested experimentally yet. Furthermore, the model postulates that the components mentioned above would play key roles. We discuss each of them in turn.

First, the models demonstrated that rhythmic patterns can be generated in a small subcircuit within each subunit of the model. The network of neurons identified extends the circuit proposed in previous work (Olivares et al., 2018), and it includes the collective participation of all motorneurons in the VNC. However, we know from experimental studies that ablation of AS motorneurons disrupts coordination of the bending wave (Tolstenkov et al., 2018). Together, this provides further evidence to the already strong hypothesis that stretch receptor feedback is a fundamental contributor to forward locomotion. It also suggests that, if the multiple intrinsic network oscillations mechanism is possible in the worm, it is a secondary, complementary contributor to its locomotion. Furthermore, if the nematode has multiple redundant mechanisms for producing locomotion (using stretch receptor feedback; using pacemaker neurons; and using network oscillators), then it follows that eliminating any one of the components in one of the three mechanisms will not cause a full disruption. To fully test the hypothesis of network rhythmic pattern generation will require analyzing the subcircuit of a VNC neural unit in isolation from proprioceptive feedback, without descending inputs from command interneurons, and in conditions where A- and B-type motorneurons are not generating intrinsic oscillations.

Second, although recent experimental work has focused on demonstrating the importance of proprioception for coordinating the multiple rhythmic pattern generators within the ventral nerve cord (Fouad et al., 2018; Gao et al., 2018; Wen et al., 2018; Xu et al., 2018), the results of our model suggest that the coordination between multiple rhythmic pattern generators can be achieved by any one of the repeating interunit electrical synapses between motorneurons. Specifically, analysis of the ensemble of models suggested three candidate electrical junctions: VB⊢⊣DB, DA⊢⊣AS, and AS⊢⊣VA. These gap junctions are present throughout the ventral nerve cord in the available connectome data. However, only a complete anatomical reconstruction of the ventral nerve cord will reveal the full extent of their presence. To test this hypothesis would require examining whether stimulation of one set of motorneurons in the ventral nerve cord can entrain the anterior or posterior motorneurons in the absence of proprioceptive feedback and with descending input from command interneurons disabled.

Finally, analysis of the ensemble of models also demonstrated that coordination between the multiple intrinsic network rhythmic pattern generators can be achieved in either the anterior or posterior direction, and some solutions coordinate in both directions. Regardless of the direction of neural pattern coordination, the behavior of the integrated neuromechanical model was a posterior traveling wave that moved forward on agar. Altogether, analysis of the solutions revealed how different configurations of the network can generate and coordinate multiple rhythmic pattern generators that will result in almost identical forward locomotion behavior.

Given the deliberate effort to remove from the model a number of components that are known to play important roles in forward locomotion for the sake of addressing a question that had not been tackled yet, the model naturally reveals discrepancies with a number of experiments. Specifically, there is evidence that: silencing A-class motor neurons does not entirely disrupt forward locomotion (Xu et al., 2018); ablating AS does not disrupt forward locomotion entirely but leads to a ventral bias in the bending (Tolstenkov et al., 2018); and D-class motor neurons may not be playing a role in coordinating body waves under certain conditions (McIntire et al., 1993; Deng et al., 2020). In contrast, the multiple network rhythmic pattern generators hypothesis in this model relies on the participation of all motor neurons, including classes A, AS, and D. Although there is further experimental work that remains to delineate with precision the role of each of those neurons in forward locomotion, there is one relatively trivial way to resolve the inconsistency between the model and the hypothesis that those neurons play no substantial roles in forward locomotion: the re-introduction of proprioceptive feedback back into the model as an independent mechanism for generating reflexive rhythmic patterns. With such a mechanism in place, silencing those neuronal classes would no longer cause disruptions to forward locomotion, as evidenced by several models that do not include those neurons, do include proprioception, and generate forward locomotion (Boyle et al., 2012; Izquierdo and Beer, 2018). That proprioceptive feedback has been removed from the model in this paper should not be taken as a postulation that it plays no role. On the contrary, we acknowledge the important role proprioceptive plays in locomotion in the worm and proceeded deliberately with the development of a model that could test the limits of whether locomotion would be possible in its absence. A second way to resolve some of the inconsistencies could arise from a complete map of the ventral nerve cord connectome. The present model only included the statistically repeating connections in the partially mapped connectome. With new connections added to the model, the possibilities for generating network rhythmic patterns only increases. With such an increase, the possibility that one of those alternatives may not require all the motor neurons in the ventral nerve cord also increases.

A more considerable discrepancy arises with the observation that body wall muscles of partially constrained worms display no activity (Supplementary Figure 1 in Wen et al., 2012). The results from those experiments suggest proprioceptive feedback is not just sufficient, but necessary for rhythmic activity in the ventral nerve cord. To the degree that proprioceptive feedback is demonstrated to be necessary for neural rhythmic activity, then the multiple network rhythmic pattern generators hypothesis must be discarded. However, despite those experiments, the possibility that the ventral nerve cord exhibits rhythmic neural patterns in the absence of proprioceptive feedback remains tenable. There is a possibility that the neural rhythmic activity generated is too weak to drive the muscles under the constrained condition. There is also a possibility that the constrained condition makes the activation of the muscles unfeasible. Although the lack of muscle activity in constrained worms provides strong evidence against intrinsic network pattern generation, it would be premature to discard the hypothesis altogether without first imaging neural activity while proprioceptive feedback is silenced.

Our model has a number of limitations. Locomotion in the worm is composed of many movements, including forward, backward, turns, and reversals (Yemini et al., 2013). A- and B- class VNC motorneurons receive descending input from command interneurons AVA and AVB, respectively, and activation from these neurons drives changes in the oscillatory patterns that lead the worm to transition between forward, backward, reversals, and turns (Chalfie et al., 1985; Kawano et al., 2011; Rakowski et al., 2013; Xu et al., 2018). Furthermore, the locomotory pattern of the worm varies drastically with different mechanical loads from the viscosity of the environment (Berri et al., 2009; Fang-Yen et al., 2010). Although a more comprehensive model would address questions regarding all of those interrelated components of locomotion in the worm, the scope of this model was much narrower. The goal of this study was not to capture all aspects of the worm's rich locomotory behavior; but to systematically search for the different parameter configurations that resulted in multiple intrinsic rhythmic pattern generators for forward locomotion on agar under the limited conditions where all VNC neurons are intact and stretch receptor feedback is not available. To study the effects of different mechanical loads on the worm's movement, it is essential to include stretch-receptor feedback into the model. To study the changes between forward and backward locomotion and changes in speed, it is essential to include the descending input from the command interneurons.

Despite the breadth of knowledge available about the worm's connectome, the signs of its connections are yet to be elucidated, and a systematic characterization of the physiological properties of the neurons and the role they play in different behaviors remains experimentally challenging, making modeling the neural basis of its behavior a massively under-constrained problem. Our approach takes this challenge seriously, generating multiple possible hypotheses for how different parameter configurations and patterns of activity could lead to the observed behavior. Despite these obvious limitations, the strength of the approach is that it integrates connectomic, physiological, and behavioral data to infer candidate configurations of synaptic properties. As experiments that map neural manipulations to behavioral kinematics increases and the data becomes public and standardized, these can be used to further constrain the optimization search and thereby hone in on the space of candidate models.

Finally, it is important to note that the fitness function used in any evolutionary experiment defines a level of constraint on the ensemble of solutions analyzed. In this work, we chose to search for model worms that could match as best as possible aspects of the worm's movement during forward locomotion on agar that have been characterized in some detail, such as the speed of movement (Cronin et al., 2005; Cohen et al., 2012). However, the conditions of the model worm are different from those of the wild type worm, crucially no stretch receptor feedback and no information from the command input neurons. The motivation of this simplification was that to the degree that a model worm without those crucial components could replicate some aspects of the wild-type behavior, then analysis of those solutions will reveal new insights about mechanisms for producing the behavior. Future work will involve setting up different levels of constraints on the model, which will result in different ensembles, each of which might be useful to understand. One interesting avenue of research can be aimed at understanding how the ensemble of solutions changes as constraints are changed.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

EO: data curation, formal analysis, supervision, investigation, visualization, methodology, and writing-original draft. EI: conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, project administration, writing-original draft, review, and editing. RB: conceptualization, resources, software, supervision, funding acquisition, validation, investigation, visualization, methodology, project administration, and writing-review and editing. All authors contributed to the article and approved the submitted version.

Funding

This material was based upon work supported by the National Science Foundation under Grant No. 1845322.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This manuscript has been released as a pre-print at https://www.biorxiv.org/content/10.1101/710566v4 (Olivares et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2021.572339/full#supplementary-material

References

Altun, Z. F., and Hall, D. H. (2009). Muscle System, Somatic Muscle. In WormAtlas. doi: 10.3908/wormatlas.1.7

CrossRef Full Text | Google Scholar

Arshavsky, Y., Deliagina, T., and Orlovsky, G. (2016). Central pattern generators: mechanisms of operation and their role in controlling automatic movements. Neurosci. Behav. Physiol. 46, 696–718. doi: 10.1007/s11055-016-0299-5

CrossRef Full Text | Google Scholar

Beer, R. D. (1995). On the dynamics of small continuous-time recurrent neural networks. Adapt. Behav. 3, 469–509. doi: 10.1177/105971239500300405

PubMed Abstract | CrossRef Full Text | Google Scholar

Berri, S., Boyle, J. H., Tassieri, M., Hope, I. A., and Cohen, N. (2009). Forward locomotion of the nematode C. elegans is achieved through modulation of a single gait. HFSP J. 3, 186–193. doi: 10.2976/1.3082260

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, J. H. (2009). C. elegans locomotion: an integrated approach (Ph.D. thesis). University of Leeds, Leeds, United Kingdom.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Butler, V. J., Branicky, R., Yemini, E., Liewald, J. F., Gottschalk, A., Kerr, R. A., et al. (2015). A consistent muscle activation strategy underlies crawling and swimming in Caenorhabditis elegans. J. R. Soc. Interface 12:20140963. doi: 10.1098/rsif.2014.0963

PubMed Abstract | CrossRef Full Text | Google Scholar

Chalasani, S. H., Chronis, N., Tsunozaki, M., Gray, J. M., Ramot, D., Goodman, M. B., et al. (2007). Dissecting a circuit for olfactory behaviour in Caenorhabditis elegans. Nature 450, 63–70. doi: 10.1038/nature06292

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiel, H. J., and Beer, R. D. (1997). The brain has a body: adaptive behavior emerges from interactions of nervous system, body and environment. Trends Neurosci. 20, 553–557. doi: 10.1016/S0166-2236(97)01149-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiel, H. J., Ting, L. H., Ekeberg, Ö., and Hartmann, M. J. Z. (2009). The brain in its body: motor control and sensing in a biomechanical context. J. Neurosci. 29, 12807–12814. doi: 10.1523/JNEUROSCI.3338-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, E., Yemini, E., Schafer, W., Feitelson, D. G., and Treinin, M. (2012). Locomotion analysis identifies roles of mechanosensory neurons in governing locomotion dynamics of C. elegans. J. Exp. Biol. 215, 3639–3648. doi: 10.1242/jeb.075416

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, N., and Boyle, J. H. (2010). Swimming at low reynolds number: a beginners guide to undulatory locomotion. Contemp. Phys. 51, 103–123. doi: 10.1080/00107510903268381

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Cohen, N., and Sanders, T. (2014). Nematode locomotion: dissecting the neuronal-environmental loop. Curr. Opin. Neurobiol. 25, 99–106. doi: 10.1016/j.conb.2013.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cook, S. J., Jarrell, T. A., Brittin, C. A., Wang, Y., Bloniarz, A. E., Yakovlev, M. A., et al. (2019). Whole-animal connectomes of both Caenorhabditis elegans sexes. Nature 571, 63–71. doi: 10.1038/s41586-019-1352-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Cronin, C. J., Mendel, J. E., Mukhtar, S., Kim, Y.-M., Stirbl, R. C., Bruck, J., et al. (2005). An automated system for measuring parameters of nematode sinusoidal movement. BMC Genetics 6:5. doi: 10.1186/1471-2156-6-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dasen, J. S. (2018). Evolution of locomotor rhythms. Trends Neurosci. 41, 648–651. doi: 10.1016/j.tins.2018.07.013

CrossRef Full Text | Google Scholar

Deng, L., Denham, J., Arya, C., Yuval, O., Cohen, N., and Haspel, G. (2020). Inhibition underlies fast undulatory locomotion in C. elegans. eNeuro. ENEURO.0241-20.2020. doi: 10.1523/ENEURO.0241-20.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, X., and Xu, J.-X. (2014). A 3D undulatory locomotion model inspired by C. elegans through DNN approach. Neurocomputing 131, 248–264. doi: 10.1016/j.neucom.2013.10.019

CrossRef Full Text | Google Scholar

Denham, J. E., Ranner, T., and Cohen, N. (2018). Signatures of proprioceptive control in Caenorhabditis elegans locomotion. Philos. Trans. R. Soc. B Biol. Sci. 373:20180208. doi: 10.1098/rstb.2018.0208

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang-Yen, C., Wyart, M., Xie, J., Kawai, R., Kodger, T., Chen, S., et al. (2010). Biomechanical analysis of gait adaptation in the nematode Caenorhabditis elegans. Proc. Natl. Acad. Sci. U.S.A. 107, 20323–20328. doi: 10.1073/pnas.1003016107

PubMed Abstract | CrossRef Full Text | Google Scholar

Faumont, S., Rondeau, G., Thiele, T. R., Lawton, K. J., McCormick, K. E., Sottile, M., et al. (2011). An image-free opto-mechanical system for creating virtual environments and imaging neuronal activity in freely moving Caenorhabditis elegans. PLoS ONE 6:e24666. doi: 10.1371/journal.pone.0024666

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenyves, B. G., Szilágyi, G. S., Vassy, Z., Sőti, C., and Csermely, P. (2020). Synaptic polarity and sign-balance prediction using gene expression data in the Caenorhabditis elegans chemical synapse neuronal connectome network. PLoS Comput. Biol. 16:e1007974. doi: 10.1371/journal.pcbi.1007974

PubMed Abstract | CrossRef Full Text | Google Scholar

Fieseler, C., Kunert-Graf, J., and Kutz, J. N. (2018). The control structure of the nematode Caenorhabditis elegans: neuro-sensory integration and proprioceptive feedback. J. Biomech. 74, 1–8. doi: 10.1016/j.jbiomech.2018.03.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Fouad, A. D., Teng, S., Mark, J. R., Liu, A., Alvarez-Illera, P., Ji, H., et al. (2018). Distributed rhythm generators underlie Caenorhabditis elegans forward locomotion. eLife 7:e29913. doi: 10.7554/eLife.29913

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, S., Guan, S. A., Fouad, A. D., Meng, J., Kawano, T., Huang, Y.-C., et al. (2018). Excitatory motor neurons are local oscillators for backward locomotion. eLife 7:e29915. doi: 10.7554/eLife.29915

PubMed Abstract | CrossRef Full Text | Google Scholar

Gjorgjieva, J., Biron, D., and Haspel, G. (2014). Neurobiology of Caenorhabditis elegans locomotion: where do we stand? Bioscience 64, 476–486. doi: 10.1093/biosci/biu058

PubMed Abstract | CrossRef Full Text | Google Scholar

Gleeson, P., Lung, D., Grosu, R., Hasani, R., and Larson, S. D. (2018). c302: a multiscale framework for modelling the nervous system of Caenorhabditis elegans. Philos. Trans. R. Soc. B Biol. Sci. 373:20170379. doi: 10.1098/rstb.2017.0379

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodman, M. B., Hall, D. H., Avery, L., and Lockery, S. R. (1998). Active currents regulate sensitivity and dynamic range in C. elegans neurons. Neuron 20, 763–772. doi: 10.1016/S0896-6273(00)81014-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Goulding, M. (2009). Circuits controlling vertebrate locomotion: moving in a new direction. Nat. Rev. Neurosci. 10:507. doi: 10.1038/nrn2608

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutenkunst, R. N., Waterfall, J. J., Casey, F. P., Brown, K. S., Myers, C. R., and Sethna, J. P. (2007). Universally sloppy parameter sensitivities in systems biology. PLoS Comput. Biol. 3:e189. doi: 10.1371/journal.pcbi.0030189

CrossRef Full Text | Google Scholar

Haspel, G., and O'Donovan, M. J. (2011). A perimotor framework reveals functional segmentation in the motoneuronal network controlling locomotion in Caenorhabditis elegans. J. Neurosci. 31, 14611–14623. doi: 10.1523/JNEUROSCI.2186-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. Proc R. Soc. Lond. B 126, 136–195. doi: 10.1098/rspb.1938.0050

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Izquierdo, E., and Beer, R. (2013). Connecting a connectome to behavior: an ensemble of neuroanatomical models of C. elegans klinotaxis. PLoS Comput. Biol. 9:e1002890. doi: 10.1371/journal.pcbi.1002890

PubMed Abstract | CrossRef Full Text | Google Scholar

Izquierdo, E., and Lockery, S. (2010). Evolution and analysis of minimal neural circuits for klinotaxis in Caenorhabditis elegans. J. Neurosci. 30, 12908–12917. doi: 10.1523/JNEUROSCI.2606-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Izquierdo, E. J. (2018). Role of simulation models in understanding the generation of behavior in C. elegans. Curr. Opin. Syst. Biol. 13, 93–101. doi: 10.1016/j.coisb.2018.11.003

CrossRef Full Text | Google Scholar

Izquierdo, E. J., and Beer, R. D. (2018). From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegans. Philos. Trans. R. Soc. B Biol. Sci. 373:20170374. doi: 10.1098/rstb.2017.0374

PubMed Abstract | CrossRef Full Text | Google Scholar

Jarrell, T. A., Wang, Y., Bloniarz, A. E., Brittin, C. A., Xu, M., Thomson, J. N., et al. (2012). The connectome of a decision-making neural network. Science 337, 437–444. doi: 10.1126/science.1221762

CrossRef Full Text | Google Scholar

Karbowski, J., Schindelman, G., Cronin, C. J., Seah, A., and Sternberg, P. W. (2008). Systems level circuit model of C. elegans undulatory locomotion: mathematical modeling and molecular genetics. J. Comput. Neurosci. 24, 253–276. doi: 10.1007/s10827-007-0054-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Katz, P. S. (2016). Evolution of central pattern generators and rhythmic behaviours. Philos. Trans. R. Soc. B Biol. Sci. 371:20150057. doi: 10.1098/rstb.2015.0057

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawano, T., Po, M. D., Gao, S., Leung, G., Ryu, W. S., and Zhen, M. (2011). An imbalancing act: gap junctions reduce the backward motor circuit activity to bias C. elegans for forward locomotion. Neuron 72, 572–586. doi: 10.1016/j.neuron.2011.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Krakauer, J. W., Ghazanfar, A. A., Gomez-Marin, A., MacIver, M. A., and Poeppel, D. (2017). Neuroscience needs behavior: correcting a reductionist bias. Neuron 93, 480–490. doi: 10.1016/j.neuron.2016.12.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Kunert, J., Proctor, J., Brunton, S., and Kutz, J. (2017). Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotion. PLoS Comput. Biol. 13:e1005303. doi: 10.1371/journal.pcbi.1005303

PubMed Abstract | CrossRef Full Text | Google Scholar

Kunert, J., Shlizerman, E., and Kutz, J. N. (2014). Low-dimensional functionality of complex network dynamics: neurosensory integration in the Caenorhabditis elegans connectome. Phys. Rev. E 89:052805. doi: 10.1103/PhysRevE.89.052805

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuramochi, M., and Doi, M. (2017). A computational model based on multi-regional calcium imaging represents the spatio-temporal dynamics in a Caenorhabditis elegans sensory neuron. PLoS ONE 12:e0168415. doi: 10.1371/journal.pone.0168415

PubMed Abstract | CrossRef Full Text | Google Scholar

Lebois, F., Sauvage, P., Py, C., Cardoso, O., Ladoux, B., Hersen, P., et al. (2012). Locomotion control of Caenorhabditis elegans through confinement. Biophys. J. 102, 2791–2798. doi: 10.1016/j.bpj.2012.04.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Leifer, A. M., Fang-Yen, C., Gershow, M., Alkema, M. J., and Samuel, A. D. (2011). Optogenetic manipulation of neural activity in freely moving Caenorhabditis elegans. Nat. Methods 8:147. doi: 10.1038/nmeth.1554

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Liu, J., Zheng, M., and Xu, X. S. (2014). Encoding of both analog-and digital-like behavioral outputs by one C. elegans interneuron. Cell 159, 751–765. doi: 10.1016/j.cell.2014.09.056

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Lindsay, T. H., Thiele, T. R., and Lockery, S. R. (2011). Optogenetic analysis of synaptic transmission in the central nervous system of the nematode Caenorhabditis elegans. Nat. Commun. 2:306. doi: 10.1038/ncomms1304

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., Hollopeter, G., and Jorgensen, E. M. (2009). Graded synaptic transmission at the Caenorhabditis elegans neuromuscular junction. Proc. Natl. Acad. Sci. U.S.A. 106, 10823–10828. doi: 10.1073/pnas.0903570106

PubMed Abstract | CrossRef Full Text | Google Scholar

Mailler, R., Avery, J., Graves, J., and Willy, N. (2010). “A biologically accurate 3d model of the locomotion of Caenorhabditis elegans,” in 2010 International Conference on Biosciences (Cancun: IEEE), 84–90. doi: 10.1109/BioSciencesWorld.2010.18

CrossRef Full Text | Google Scholar

Marder, E., and Bucher, D. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, R986-R996. doi: 10.1016/S0960-9822(01)00581-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Marom, S., and Shahaf, G. (2002). Development, learning and memory in large random networks of cortical neurons: lessons beyond anatomy. Q. Rev. Biophys. 35:63. doi: 10.1017/S0033583501003742

PubMed Abstract | CrossRef Full Text | Google Scholar

McIntire, S. L., Jorgensen, E., Kaplan, J., and Horvitz, H. R. (1993). The GABAergic nervous system of Caenorhabditis elegans. Nature 364:337. doi: 10.1038/364337a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellem, J. E., Brockie, P. J., Madsen, D. M., and Maricq, A. V. (2008). Action potentials contribute to neuronal signaling in C. elegans. Nat. Neurosci. 11, 865–867. doi: 10.1038/nn.2131

PubMed Abstract | CrossRef Full Text | Google Scholar

Milligan, B., Curtin, N., and Bone, Q. (1997). Contractile properties of obliquely striated muscle from the mantle of squid (Alloteuthis subulata) and cuttlefish (Sepia officinalis). Journal of Exp. Biol. 200, 2425–2436.

PubMed Abstract | Google Scholar

Minassian, K., Hofstoetter, U. S., Dzeladini, F., Guertin, P. A., and Ijspeert, A. (2017). The human central pattern generator for locomotion: does it exist and contribute to walking? Neuroscientist 23, 649–663. doi: 10.1177/1073858417699790

PubMed Abstract | CrossRef Full Text | Google Scholar

Niebur, E., and Erdös, P. (1991). Theory of the locomotion of nematodes: dynamics of undulatory progression on a surface. Biophys. J. 60, 1132–1146. doi: 10.1016/S0006-3495(91)82149-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Olivares, E., Izquierdo, E. J., and Beer, R. D. (2020). A neuromechanical model of multiple network rhythmic pattern generators for forward locomotion in C. elegans. bioRxiv. doi: 10.1101/710566

CrossRef Full Text | Google Scholar

Olivares, E. O., Izquierdo, E. J., and Beer, R. D. (2018). Potential role of a ventral nerve cord central pattern generator in forward and backward locomotion in Caenorhabditis elegans. Netw. Neurosci. 2, 323–343. doi: 10.1162/netn_a_00036

PubMed Abstract | CrossRef Full Text | Google Scholar

Palyanov, A., Khayrulin, S., and Larson, S. D. (2018). Three-dimensional simulation of the Caenorhabditis elegans body and muscle cells in liquid and gel environments for behavioural analysis. Philos. Trans. R. Soc. B Biol. Sci. 373:20170376. doi: 10.1098/rstb.2017.0376

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastore, V. P., Massobrio, P., Godjoski, A., and Martinoia, S. (2018). Identification of excitatory-inhibitory links and network topology in large-scale neuronal assemblies from multi-electrode recordings. PLoS Comput. Biol. 14:e1006381. doi: 10.1371/journal.pcbi.1006381

PubMed Abstract | CrossRef Full Text | Google Scholar

Peliti, M., Chuang, J. S., and Shaham, S. (2013). Directional locomotion of C. elegans in the absence of external stimuli. PLoS ONE 8:e78535. doi: 10.1371/journal.pone.0078535

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfeifer, R., and Bongard, J. (2006). How the Body Shapes the Way We Think: A New View of Intelligence. MIT Press. doi: 10.7551/mitpress/3585.001.0001

CrossRef Full Text | Google Scholar

Pierce-Shimomura, J. T., Chen, B. L., Mun, J. J., Ho, R., Sarkis, R., and McIntire, S. L. (2008). Genetic analysis of crawling and swimming locomotory patterns in C. elegans. Proc. Natl. Acad. Sci. U.S.A. 105, 20982–20987. doi: 10.1073/pnas.0810359105

PubMed Abstract | CrossRef Full Text | Google Scholar

Prinz, A., Bucher, D., and Marder, E. (2004). Similar network activity from disparate circuit parameters. Nat. Neurosci. 7, 1345–1352. doi: 10.1038/nn1352

PubMed Abstract | CrossRef Full Text | Google Scholar

Putrenko, I., Zakikhani, M., and Dent, J. A. (2005). A family of acetylcholine-gated chloride channel subunits in Caenorhabditis elegans. J. Biol. Chem. 280, 6392–6398. doi: 10.1074/jbc.M412644200

PubMed Abstract | CrossRef Full Text | Google Scholar

Rakowski, F., Srinivasan, J., Sternberg, P., and Karbowski, J. (2013). Synaptic polarity of the interneuron circuit controlling C. elegans locomotion. Front. Comput. Neurosci. 7:128. doi: 10.3389/fncom.2013.00128

PubMed Abstract | CrossRef Full Text | Google Scholar

Rand, J. B., Duerr, J. S., and Frisby, D. L. (1997). “Neurogenetics of synaptic transmission in Caenorhabditis elegans,” in Advances in Pharmacology, Vol. 42, eds D. S. Goldstein, G. Eisenhofer, and R. McCarty (Academic Press), 940–944. doi: 10.1016/S1054-3589(08)60902-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossignol, S., Dubuc, R., and Gossard, J.-P. (2006). Dynamic sensorimotor interactions in locomotion. Physiol. Rev. 86, 89–154. doi: 10.1152/physrev.00028.2005

CrossRef Full Text | Google Scholar

Sauvage, P. (2007). Etude de la locomotion chez C. elegans et perturbations mecaniques du mouvement (Ph.D. thesis). Laboratoire Matiere et Systemes Compiexes, Paris, France.

Google Scholar

Tavernarakis, N., Shreffler, W., Wang, S., and Driscoll, M. (1997). unc-8, a DEG/ENaC family member, encodes a subunit of a candidate mechanically gated channel that modulates C. elegans locomotion. Neuron 18, 107–119. doi: 10.1016/S0896-6273(01)80050-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Tolstenkov, O., Van der Auwera, P., Costa, W. S., Bazhanova, O., Gemeinhardt, T. M., Bergs, A. C., et al. (2018). Functionally asymmetric motor neurons contribute to coordinating locomotion of Caenorhabditis elegans. eLife 7:e34997. doi: 10.7554/eLife.34997

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshney, L., Chen, B., Paniagua, E., Hall, D., and Chklovskii, D. (2011). Structural properties of the Caenorhabditis elegans neuronal network. PLoS Comput. Biol. 7:e1001066. doi: 10.1371/journal.pcbi.1001066

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidal-Gadea, A., Topper, S., Young, L., Crisp, A., Kressin, L., Elbel, E., et al. (2011). Caenorhabditis elegans selects distinct crawling and swimming gaits via dopamine and serotonin. Proc. Natl. Acad. Sci. U.S.A. 108, 17504–17509. doi: 10.1073/pnas.1108673108

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, H. (1969). Wave formation by infective larvae of the plant parasitic nematode meloidogyne javanica. Nematologica 15, 65–75. doi: 10.1163/187529269X00100

CrossRef Full Text | Google Scholar

Waterston, R. (1988). “Muscle,” in The Nematode C. elegans, ed W. Wood (New York, NY: Cold Spring Harbor Laboratory Press), 281–335.

Google Scholar

Wen, Q., Gao, S., and Zhen, M. (2018). Caenorhabditis elegans excitatory ventral cord motor neurons derive rhythm for body undulation. Philos. Trans. R. Soc. B Biol. Sci. 373:20170370. doi: 10.1098/rstb.2017.0370

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, Q., Po, M. D., Hulme, E., Chen, S., Liu, X., Kwok, S. W., et al. (2012). Proprioceptive coupling within motor neurons drives C. elegans forward locomotion. Neuron 76, 750–761. doi: 10.1016/j.neuron.2012.08.039

PubMed Abstract | CrossRef Full Text | Google Scholar

White, J., Southgate, E., Thomson, J., and Brenner, S. (1986). The structure of the nervous system of the nematode Caenorhabditis elegans. Philos. Trans. R. Soc B Biol. Sci. 275, 327–348.

PubMed Abstract | Google Scholar

Wicks, S. R., Roehrig, C. J., and Rankin, C. H. (1996). A dynamic network simulation of the nematode tap withdrawal circuit: predictions concerning synaptic function using behavioral criteria. J. Neurosci. 16, 4017–4031. doi: 10.1523/JNEUROSCI.16-12-04017.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, M., Jarrell, T. A., Wang, Y., Cook, S. J., Hall, D. H., and Emmons, S. W. (2013). Computer assisted assembly of connectomes from electron micrographs: application to Caenorhabditis elegans. PLoS ONE 8:e54050. doi: 10.1371/journal.pone.0054050

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, T., Huo, J., Shao, S., Po, M., Kawano, T., Lu, Y., et al. (2018). Descending pathway facilitates undulatory wave propagation in Caenorhabditis elegans through gap junctions. Proc. Natl. Acad. Sci. U.S.A. 115, E4493-E4502. doi: 10.1073/pnas.1717022115

PubMed Abstract | CrossRef Full Text | Google Scholar

Yemini, E., Jucikas, T., Grundy, L. J., Brown, A. E., and Schafer, W. R. (2013). A database of Caenorhabditis elegans behavioral phenotypes. Nat. Methods 10:877. doi: 10.1038/nmeth.2560

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeon, J., Kim, J., Kim, D.-Y., Kim, H., Kim, J., Du, E. J., et al. (2018). A sensory-motor neuron type mediates proprioceptive coordination of steering in C. elegans via two trpc channels. PLoS Biol. 16:e2004929. doi: 10.1371/journal.pbio.2004929

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhen, M., and Samuel, A. D. (2015). C. elegans locomotion: small circuits, complex functions. Curr. Opin. Neurobiol. 33, 117–126. doi: 10.1016/j.conb.2015.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar