GitHub pages
The nematode Caenorhabditis elegans possesses a remarkably small but fully mapped nervous system, consisting of precisely 302 neurons. This simplicity, combined with the availability of its complete connectome, makes it the most tractable multicellular organism for computational modeling of neural dynamics. The first comprehensive connectome of C. elegans was reconstructed by White and colleagues in 1986 using serial electron microscopy, representing a milestone in neuroscience. This map revealed the identities of individual neurons, their types, and the complete set of chemical and electrical synaptic connections between them.
The neurons in the connectome are labeled using a standardized naming convention that encodes anatomical position and type. For instance, a neuron labeled ADEL indicates an anterior dorsal left neuron, while ADAL indicates an anterior dorsal left neuron of a slightly different type. These names serve as unique identifiers in the connectome dataset. Importantly, neurons may form multiple types of connections. Chemical synapses are directional and typically represented as an asymmetric adjacency matrix, whereas electrical synapses (gap junctions) are bidirectional and represented by a symmetric adjacency matrix.
To mathematically model the network dynamics, we define two adjacency matrices: $A_{\text{chem}}$ for chemical synapses and $A_{\text{gap}}$ for gap junctions. These matrices may be combined linearly to form a single weighted adjacency matrix:
\[A_{\text{comb}} = A_{\text{chem}} + k \cdot A_{\text{gap}},\]where $k$ is a scalar coefficient that adjusts the relative contribution of gap junctions. Each element $A_{\text{comb}}[j,i]$ represents the strength of the connection from neuron $i$ (presynaptic) to neuron $j$ (postsynaptic). The rows correspond to postsynaptic neurons, and the columns correspond to presynaptic neurons. If neuron $i$ connects to multiple downstream neurons, each postsynaptic neuron receives an input scaled by the respective synaptic weight $A[j,i]$.
A crucial step in analyzing the connectome is identifying candidate sensory and motor neurons. One commonly used heuristic is based on the network degrees. The outdegree of neuron $i$ is defined as the sum of its outgoing connections:
\[\text{outdeg}[i] = \sum_{j=1}^{N} A_{\text{comb}}[j,i],\]while the indegree is the sum of its incoming connections:
\[\text{indeg}[i] = \sum_{j=1}^{N} A_{\text{comb}}[i,j].\]The difference between outdegree and indegree, $\text{diff}[i] = \text{outdeg}[i] - \text{indeg}[i]$, provides a heuristic for functional role. Neurons with large positive differences are likely to be sensory neurons, as they predominantly send information to the network, whereas neurons with large negative differences are likely to be motor neurons, as they predominantly receive information from the network.
The dynamics of each neuron in a rate-based model are represented by a continuous variable $x_i(t)$, which abstracts the activity of neuron $i$. This activity can be interpreted as the instantaneous firing rate or the membrane potential of the neuron. The rate-based neuron obeys the differential equation:
\[\frac{dx_i}{dt} = -\frac{x_i}{\tau} + \sum_{j=1}^{N} A_{ij} f(x_j) + I_i,\]where $N$ is the total number of neurons, $A_{ij}$ is the element of the adjacency matrix corresponding to the connection from neuron $j$ to neuron $i$, $f(x_j)$ is a possibly nonlinear activation function applied to presynaptic neuron $j$, $I_i$ is an externally applied input to neuron $i$, and $\tau$ is the membrane time constant.
Each term in the model equation has a direct biological interpretation.
The variable $x_i(t)$ represents the internal state or activity of neuron $i$. In biological terms, this could correspond to the membrane potential or the firing rate of the neuron.
The first term, $-x_i/\tau$, represents the exponential decay of the neuron’s activity toward a baseline, modeling the passive leakage of membrane potential across the neuronal membrane.
The membrane time constant $\tau$ is given by the product of the membrane resistance $R_m$ and the membrane capacitance $C_m$, determining how quickly the neuron responds to input and how long it retains past activity.
A small $\tau$ produces rapid decay and fast responses, whereas a large $\tau$ produces slower decay and longer temporal integration.
The second term, $\sum_{j} A_{ij} f(x_j)$, represents the total synaptic input to neuron $i$ from all presynaptic neurons. Each presynaptic neuron $j$ contributes $f(x_j) \cdot A_{ij}$ to the postsynaptic input, where $f(x_j)$ is the output of neuron $j$ after applying an activation function. If the activation function is nonlinear, such as $\tanh$, it simulates biological saturation: the output is bounded, reflecting limits on firing rate or postsynaptic potential. If the activation function is the identity, the model reduces to a linear system. Gap junctions can be included in $A_{ij}$ with an appropriate scaling factor to represent the electrical coupling between neurons.
The third term, $I_i$, represents any external input applied to the neuron, which can model sensory stimuli in biological terms or experimental current injection. This allows selective activation of sensory neurons to probe network response.
The presynaptic contribution is weighted, meaning that not all postsynaptic neurons receive the same input from a firing neuron. Strong synapses produce a larger influence on the postsynaptic neuron, analogous to greater neurotransmitter release or higher conductance at a gap junction in biological neurons. The output signal of a neuron is thus a combination of its internal state, synaptic weights, and the nonlinear activation function, and it represents an abstract “strength of influence” rather than a literal spike count.
The equation is integrated numerically using the Euler method. Given a small timestep $\Delta t$, the update rule is:
\[x_i(t + \Delta t) = x_i(t) + \Delta t \left( -\frac{x_i(t)}{\tau} + \sum_{j} A_{ij} f(x_j(t)) + I_i(t) \right).\]This produces a time series of neuron activities that can be analyzed to study information propagation, network dynamics, and the response of motor neurons to sensory stimuli.
Normalization is often applied to the adjacency matrix to prevent unbounded growth of activity:
\[\tilde{A}_{ij} = \frac{A_{ij}}{\sum_{k} |A_{kj}| + \epsilon},\]where $\epsilon$ is a small number added to prevent division by zero for neurons with no outgoing connections.
Historically, rate-based models have been used since the 1990s as a tractable approximation of neural dynamics before moving to spiking neuron models. They abstract away individual action potentials while retaining the essence of synaptic interactions and network effects, allowing simulations of full connectomes like that of C. elegans. Nonlinear activation functions allow the model to capture saturation and thresholding effects observed in biological neurons, while the decay term represents passive leakage of membrane potential, giving neurons an intrinsic timescale of response.
In conclusion, modeling the C. elegans connectome using a rate-based framework involves representing neuron activity as continuous variables, weighted synaptic interactions, decay toward baseline activity, and optional nonlinear activation. The model abstracts the essential features of neural computation while remaining computationally tractable, providing a bridge between biological neuron properties and mathematical network analysis.
The C. elegans connectome is typically stored as a tabular dataset, where each row corresponds to a connection between a presynaptic neuron and a postsynaptic neuron. A simplified excerpt of such a dataset might look like:
| Presynaptic Neuron | Postsynaptic Neuron | Connection Type | Count |
|---|---|---|---|
| ADAR | ADAL | EJ (Electrical Junction) | 1 |
| ADFL | ADAL | EJ | 1 |
| ASHL | ADAL | EJ | 1 |
| AVDR | ADAL | EJ | 2 |
| PVQL | ADAL | EJ | 1 |
| ADEL | ADAL | Sp (Chemical Synapse) | 1 |
| ADFL | ADAL | Sp | 1 |
| AIAL | ADAL | Sp | 1 |
| AIBL | ADAL | R (Reciprocal Synapse) | 1 |
| AIBR | ADAL | Rp (Reciprocal) | 2 |
A = Anterior (toward the head)P = Posterior (toward the tail)L = LeftR = RightSp = chemical synapse (directional)EJ = gap junction (bidirectional)R / Rp = reciprocal connectionsFrom this raw data, the adjacency matrix is constructed as follows:
Chemical synapse matrix, $A_{\text{chem}}$: Each entry $A_{\text{chem}}[j,i]$ is the sum of all chemical synapses from neuron $i$ to neuron $j$. The matrix is generally asymmetric.
Gap junction matrix, $A_{\text{gap}}$: Each entry $A_{\text{gap}}[j,i]$ counts the electrical junctions. The matrix is symmetric because current flows in both directions.
After optionally scaling gap junctions, the combined adjacency matrix used in the rate-based model is:
\[A_{\text{comb}} = A_{\text{chem}} + k \cdot A_{\text{gap}}\]Below is a schematic of a neuron in the network, highlighting all terms from the rate-based equation and their biological counterparts:
External Input I_i
│
▼
┌───────────┐
Presynaptic │ │ Postsynaptic
│ Neuron i │ ┌─────────┐ ┌─────────┐
f(x_1)───┐ │ x_i │──▶ │ Neuron j│ │ Neuron k│
f(x_2)───┼──▶│ │──▶ │ x_j │ │ x_k │
f(x_3)───┘ │ │ └─────────┘ └─────────┘
└───────────┘
▲
│
Decay Term
-x_i/τ
Presynaptic Inputs: The neuron receives contributions from other neurons, scaled by the adjacency matrix A[j,i] and optionally passed through a nonlinear function f(x_j). Different postsynaptic neurons receive different inputs depending on synaptic weights.
Decay Term: The neuron’s activity naturally decays toward baseline at a rate determined by the membrane time constant $\tau$.
External Input: Sensory neurons or experimentally stimulated neurons receive an additional input $I_i$.
Postsynaptic Outputs: The activity of neuron i, after optional nonlinear transformation, drives downstream neurons. This represents the influence that neuron i exerts on its targets.
This schematic illustrates how each component of the rate-based differential equation maps directly onto network connectivity and biological processes.
Here’s my rudimentay implementation of the model in Python using PyTorch: https://github.com/pankajgarkoti/elegans-neural-network
However, there are a few established projects that aim to completely simulate the entire working of this creature from the connectome data. Please refer to the following for a complete implementation: