A minimal reaction-diffusion neural model generates C. elegans undulation

The small (1 mm) nematode Caenorhabditis elegans (C. elegans) has become a widely used model organism (cf. http://www.wormbook.org [1]), and has been among the most studied biological models of neuronal development and locomotion [5, 1]. The C. elegans connectome has been completely mapped [6] and, as described below, its locomotion has been widely studied (c.f. [1, 7]). There are a variety of neuronal models which can generate such undulation, “When crawling on a solid surface, the nematode C. elegans moves forward by propagating sinusoidal dorso-ventral retrograde contraction waves. A uniform propagating wave leads to motion that undulates about a straight line.” [8]. A different type of locomotion, often called swimming, occurs when nematodes are submerged in a liquid medium. The nematodes “switch” between these two gaits, by changing the dynamics of the central pattern generator (CPG).

The purpose of this paper is to describe a minimal, biomimetic, reaction-diffusion model for the C. elegans central pattern generator (CPG) [2, 9]. We use simulation methods to show that a small network of [10]-[11] neurons, see also [12] (one of the simplest neuronal models), based on a skeleton model of the C. elegans CPG, can reproduce key features of C. elegans undulation [2] [13] [14], and thus locomotion.

Finally, we describe an analog electronic implementation of our model through solving a modified version of the FitzHugh-Nagumo neuron [10], based upon an analog circuit originally proposed by [3]. This circuit solves the Keener differential equations, and we adjusted it to allow diffusive coupling between neurons. We constructed a small network with these “neuro-mimetic” circuits, and showed that their behaviour replicates FitzHugh-Nagumo simulated behaviour.

There are many other CPG models; for example, [4] proposes a distributed network of self-oscillating systems of neurons, instead of a structured chain like our proposed CPG. Here we describe a minimal working model, rather than striving for fuller realism, in order to explore the fundamental components of a small CPG. We aimed for a simple, ”minimal” model to enable exploration of defects or other changes in the CPG in an efficient, reproducible and explainable way. Our model can thus be considered as intermediate between the early ”neuro-mechanical” model of [15] and more realistic models such as [4], a role similar to that of [16] biochemical computation.

Recall that central pattern generator is a small neural circuit which generates and regulates the movement of complex organisms. This structure is present in different forms in many animals, and regulates many types of periodic motion. Changes in gait are driven by changes in the dynamics of the CPG, c.f. [17]. The CPG circuit topology is constant; the mode of locomotion depends on the sequence in which the neurons fire.

Our simple model C. elegans central pattern generator has two principal components. The first is the head oscillator. As described by [7], the head oscillator consists of two “head neurons” with mutually inhibitory coupling. Oscillations are generated when this coupling destabilizes an excitable steady state. Here two Fitzhugh-Nagumo neurons, with oscillatory dynamics, stabilized $180^{\circ}$ out of phase by mutually inhibitory coupling. This provides a pair of out-of-phase stimuli which propagate through a descending pathway of pairs of coupled, excitable, dorsal and ventral neurons. These follow the body of the worm, and are linked to motor neurons and muscles. The head oscillator drives the descending pathway, and the pathway is kept in sync by mutual inhibitory coupling between neurons.

We use 12 pairs of such neurons, as in C. elegans coupled by model gap junctions. This coupling yields phase lags as we descend the pathway from head to tail, enabling the head oscillator to drive traveling waves of excitation along the body of the nematode. See \Freffig: xu_cpg for descending pathway of [2] and our simplified model; the latter depicted as a graph, wherein neurons are nodes, and the arrows between them represent connections.

As described above, we sought to use the simplest relevant neuronal model. The classical Hodgkin-Huxley[21] model of squid neurons has led to a variety of simpler conduction models, including the Morris-Lecar[22] and [10]-[11] (FHN) models. The FHN model consists of two dynamical variables; a fast activator variable $v$ corresponding to the (rescaled) membrane potential, and a slow inhibitor variable $w$ corresponding to a generalized gating variable.

The parameter $I_{\mathrm{ext}}$ is an external driving current, and is used here to model the effect of gap junctions and synapses upon membrane potential. The parameter $\alpha$ determines the vertical position of the $w$-nullcline and thus controls excitability. Action potentials can also be generated by a current injection corresponding to $I_{\mathrm{ext}}$. Finally, the parameter $\varepsilon$ determines ratio between the time scales of the fast activator variable $v$ and the slow inhibitor variable $w$; $\gamma$ determines the effect of the growth of the slow gate variable $w$. We used the following standard parameter values: $\varepsilon=0.08,\gamma=0.8$ [12]. FitzHUgh-Nagumo neurons can display either (self-)oscillatory or excitable dynamics; given a sufficiently large stimulus (here a pulse of injected current $I$), an excitable neuron generates an action potential. An oscillatory neuron generates a regular sequence of action potentials. In comparison to the standard value $\alpha=0.7$ [12], we adjusted the parameter $\alpha$ to control neuronal dynamics, with values larger than the oscillatory-excitable boundary $\alpha_{0}\approx 0.467$ generating excitable dynamics from a stable steady state, and smaller values generating oscillatory dynamics as the steady state is destabilized.

Moreover, $f(v)$ can be any function which retains the appropriate dynamics, in that it has the same general shape as the cubic $f(v)=\frac{v^{3}}{3}-v$. For example, $f(v)$ could be replaced by the cubic-like I-V curve of the tunnel diode in the [11] circuit, or even the piecewise linear approximation generated in the [3] circuit used in our implementation.

[2] described a simplified two-variable model for C. elegans neurons, consisting of a fast, cubic-like activator variable (see the $v$-nullcline in \Freffig: nullclines) and a slow, non-linear inhibitor variable (see the $n$-nullcline). The [2] neuronal model can thus be interpreted as a Fitzhugh-Nagumo type model neuron.

Following [17, 23], we now use generalized diffusion coupling between FHN neurons to model gap junctions and synapses. The effect of gap junctions and synapses upon the membrane potential is modeled by replacing the external current $I_{\mathrm{ext}}$ of the system, by a generalized diffusion term, $D\max(\Delta v,0)$. A positive diffusion coefficient $D$ is used to simulate a gap junction, or an excitatory synapse and a negative coefficient to simulate an inhibitory synapse ([17]).

This yields a system of diffusion-coupled FitzHugh-Nagumo equations, c.f. [24]:

where $\Delta v$ is the rectified difference in voltage between the driving and driven neurons, essentially $\Delta v=max(v_{driven}-v_{driving},0)$. For example, consider the potential of an oscillatory FitzHugh-Nagumo neuron to drive an excitable FitzHugh-Nagumo neuron. The consequent dynamics depends upon the magnitude of positive diffusion coupling $D>0$: $D$ must be sufficiently large to generate an action potential in the driven excitable neuron, which then responds with a time delay correspodning inversely to the magnitude of $D$.

We start with the network shown in \Freffig: xu_cpg, bottom. We perform simulations in Python, using the standard SciPy ODE solvers, which wrap LSODA, the Livermore Solver for Ordinary Differential Equations. The following block of code illustrates these calculations. Much like [25] did, we take a segment-based approach to the worm model. Pairs of muscles on either side determine the overall angular displacement per segment. The magnitude of the “contraction” resulting from the smoothing is interpreted as an angular displacement for that segment. We use cubic spline interpolation to smooth the worm body.

Muscle response is simulated by applying Gaussian filters to the positive component of membrane potentials:

The worm is then described initially as a sequence of nodes, joined by line segments, with the angle at each node determined by muscular tension. This is followed by a cubic spline interpolation. Finally, we generated a simple video by fixing the head of a worm to the origin; see \Freffig: worm_neuron_dash.

Simulation versus experimental results. Our model generates traveling waves of approximately sinusoidal form, with approximately constant amplitude from head to tail. In comparison, real C. elegans undulations are characterized by approximately sinusoidal oscillations, but with amplitude decreasing slightly from head to tail, c.f. [2]. In addition, a search for ”markers of chaos” as in [14] finds that undulations of the simulated worm are too regular (neutrally stable), in contrast to unduations of real nematodes (with $\mathrm{Max}\left\{Re(\lambda)\right\}\approx 1\mathrm{s}^{-1}$), likely because we have yet to implement proprioceptive neuronal inputs, c.f. [26] and references therein. The neutral stability of our model may allow it to respond readily but relatively consistently to inputs. Moreover, a small, negative $Max\{Re(\lambda)\}$ would typically maintain relatively consistent undulation, but be readily overcome by somewhat larger inputs or noise, the latter generating random changes in undulation pattern.

Further comparisons and a more detailed non-linear analysis of real C. elegans undulation will appear in a forthcoming paper [Zhang et al., in preparation].

We show that a diffusion-coupled [3] electronic analog neurons can imitate key features of the dynamics of our network of Fitzhugh-Nagumo neurons. Recall that [11] proposed a circuit to simulate a FitzHugh-Nagumo neuron, shown in \Freffig: analog_neurons. It used a tunnel diode to achieve a cubic-like activation function, and an inductor to differentiate the potassium current (slow gate variable) \Freffig: analog_neurons

However, given that tunnel diodes are expensive and rarely available, [3] proposed a modified Nagumo circuit which used the saturation properties of operational amplifiers (”op-amps”) to achieve cubic-like non-linearity in the FHN model. This yields a piecewise linear approximation to the Fitzhugh-Nagumo $v$ (sodium, fast activator) nullclines, as shown in \Freffig: nullclines. The nullclines are sufficiently similar that the dynamics are effectively the same [3]. Keener also used an op-amp and a capacitor to simulate the inductor in the original Nagumo circuit. Roughly, differentiation is simulated by ”anti-integration.Later [27] proposed a CMOS implementation of the Nagumo circuit.

Finally, diffusion coupling is simulated by a follower, and then a coupling resistor followed by a diode (positive diffusion, excitatory synapse or gap junction) or a follower, inverting amplifier of unit gain, coupling resistor and a diode in series (negative diffusion, inhibitory synapse), c.f. [17, 23] for the role of diffusion, c.f. [24] for the use of resistors to emulate diffusion. The diffusion constant is inversely proportional to the coupling resistance. Coupling resistors are tuned to obtain desired synchronization or time delays; for example, coupling resistors in the descending chain of neurons are $\sim$ 1-3 M$\Omega$, with larger resistors yielding longer delays until eventually the injected current is too small to generate an action potential. Typical experimental results for a pair of neurons in the descending chain coupled by a gap junction are shown in \Freffig: Keener_simulation.

We have shown that the undulatory motion of C. elegans can be generated using a simple network of Fitzhugh-Nagumo neurons, the network capturing using a structured central pattern generator, and simple, biomimetic neurons. In particular, our simple, ”minimal” model generates traveling waves of sinusoidal undulations, as in [2]. An initial exploration shows that small changes in parameters do not alter qualitative dynamics; further study is needed to explore the range of possible dynamics.

Our model thus lies in between the ”Shape memory alloy-based small crawling robot” mechanical model of [15] and more realistic (and complex) descriptions of [25] and [4], in the spirit of [16] description of the neural system as a computational system. Because of the simplicity and flexibility our our model, we expect that it may prove useful is studying the effects of stimuli, aging and mutations upon the dynamics of C. elegans undulation and locomotion.

Acknowledgement. We acknowledge assistance from our students Cheris Congo, Miranda Hulsey-Vincent, Rifah Tasnim and Naol Negassa.