Searching for collective behavior in a small brain

In large neuronal networks, it is believed that functions emerge through the collective behavior of many interconnected neurons. Recently, the development of experimental techniques that allow simultaneous recording of calcium concentration from a large fraction of all neurons in Caenorhabditis elegans—a nematode with 302 neurons—creates the opportunity to ask whether such emergence is universal, reaching down to even the smallest brains. Here, we measure the activity of 50+ neurons in C. elegans, and analyze the data by building the maximum entropy model that matches the mean activity and pairwise correlations among these neurons. To capture the graded nature of the cells' responses, we assign each cell multiple states. These models, which are equivalent to a family of Potts glasses, successfully predict higher statistical structure in the network. In addition, these models exhibit signatures of collective behavior: the state of single cells can be predicted from the state of the rest of the network; the network, despite being sparse in a way similar to the structural connectome, distributes its response globally when locally perturbed; the distribution over network states has multiple local maxima, as in models of memory; and the parameters that describe the real network are close to a critical surface in this family of models.

Figure 1

Schematics of data acquisition and processing. (a) Examples of the raw images acquired through the $10\times$ (scale bar equals $100\mu m$) and $40\times$ (scale bar equals $10\mu m$) objectives. The body of the nematode is outlined with light green curves. As an example, we show for one neuron that (b) the intensity of the nuclei-localized fluorescent protein tags—the calcium-sensitive GCaMP and the control fluorophore RFP—are measured as functions of time. Photobleaching occurs on a longer timescale than the intracellular calcium dynamics, which allows us to perform photobleaching correction by dividing the raw signal with its exponential fit, resulting in the signals of panel (c). (d) The normalized ratio of the photobleaching-corrected intensity, $f$, is a proxy for the calcium concentration in each neuron nuclei (dark gray). As described in the text, this signal is discretized using the denoised time derivative $\dot{f}$; we use three states, marked as red, blue, and black after smoothing (lightly offset for ease of visualization). (e) The time derivative $\dot{f}$, extracted using total-variation regularized differentiation.

Figure 2

Comparison of pairwise mutual information distribution for the calcium-sensitive GCaMP worms and the GFP control worms. Mutual information is estimated using binning and finite-sample extrapolation methods as described in [35] for all pairs of neurons. For the normalized fluorescence ratio, $f$, the distribution of the mutual information, $P(MI(f_i;f_j))$, exhibits little difference between the calcium-sensitive GCaMP worm and the GFP control worm [panel (a)]. In comparison, for the time derivative of the normalized fluorescence ratio, $\dot{f}$, the distribution of the mutual information, $P(MI(\dot{f_i};\dot{f_j}))$, is peaked around zero for the GFP control worm, while the distribution is wide for the calcium-sensitive GCaMP worm [panel (b)]. This observation suggests that the time derivative of the fluorescence ratio, $\dot{f_i}$, is more informative than its magnitude, $f_i$.

Figure 3

Discretization of the empirically observed fluorescence signals. (a) Heat map of the normalized fluorescence ratio between photobleaching-corrected GCaMP fluorescence intensity and RFP fluorescence intensity, $f$, for each neuron as a function of time. (b) Heat map of the neuronal activity after discretization based on time derivatives of $f$. Green corresponds to a state of “rising,” red “falling,” and white “flat.”

Figure 4

Model construction: Learning the maximum entropy model from data. (a) Connected pairwise correlation matrix, $C_{ij}$, measured for a subgroup of 50 neurons. (b) The inferred interaction matrix, $J_{ij}$. (c) Probability of neuron $i$ in state $r$ (green circles for the “rise” state, red crosses for the “fall” state, and black dots for the “flat” state), for the same group of 50 neurons as panel (a). (d) The inferred local field, $h_i^r$. (e) Model reproduces pairwise correlation (unconnected) within variation throughout the experiment. Error bars are extrapolated from bootstrapping random halves of the data. (f) Same as panel (e), but for mean neuron activity $m_i^r$.

Figure 5

(a) No signs of overfitting are observed for pairwise maximum entropy models with up to $N=50$ neurons, measured by the difference of per-neuron log likelihood of the data under the pairwise maximum entropy model for training sets consisting of $5/6$ of the data and test sets. Clusters around $N=10,15,20,\cdots ,50$ represent randomly chosen subgroups of $N$ neurons. Error bars are the standard deviation across 10 random partitions of training and test samples. The dashed lines show the expected per-neuron log-likelihood difference and its standard deviation calculated through perturbation methods (see Appendix pp1). (b) The difference between log likelihood of the training data and of the test data is greater than 0 (the red line) within error bars for maximum entropy models on $N=10,20,\cdots ,50$ neurons with pairwise correlation tensor constraint (see Appendix pp2), which suggests that this model does not generalize well.

Figure 6

Model validation: The model predicts unconstrained higher order correlations of the data. Panel (a) shows the comparison between model prediction and data for the connected three-point correlation $C_{ijk}$ for a representative group of $N=50$ neurons. All 19 800 possible triplets are plotted with the blue dot. Error bars are generated by bootstrapping random halves of the data, and are shown for 20 uniformly spaced random triplets in red. Panel (b) shows the error of three-point function $\mathrm{\Delta }{C}_{ijk}$ as a function of the connected three-point function $C_{ijk}$, binned by its value predicted by the $C_{ijk}$ model. The red curve is the difference between data and model prediction. The blue curve is the standard error from mean of $C_{ijk}$ over the course of the experiment, extracted by bootstrapping random halves of the experiment. Panels (c) and (d) are the same as panels (a) and (b), but for the connected two-point correlation tensor $C_{ij}^{rs}$.

Figure 7

Model validation: Comparison between model prediction and data for observables not constrained by the model. The neuron network has $N=50$ neurons. (a) Probability of $k$ neurons being in the same state. Blue dots are computed from the data. Yellow dash-dotted line is the prediction from a model in which all neurons are independent, generated by applying a random temporal cyclic permutation to the activity of each neuron. Purple line is the prediction of the pairwise maximum entropy model. (b) Tail distribution of the energy for the data and the model. All error bars in this figure are extrapolated from bootstrapping. (c), (d) Probability ratio of the state of a single neuron as a function of the effective field $g_i^r$, binned by the value of the effective field. Error bars are the standard deviation after binning.

Figure 8

Energy landscape of the inferred maximum entropy model. (a) Schematic of the energy landscape with local minima $\alpha$, $\beta$ and the corresponding basin ${\mathrm{\Omega }}_{\alpha }$, ${\mathrm{\Omega }}_{\beta }$. Colored in light blue is the metabasin formed at the given energy threshold, $\mathrm{\Delta }E$. (b) Typical distribution of the value of the energy minima and the barriers of a maximum entropy model on $N=30$ neurons. The global energy minimum, $E_0$, is subtracted from the energy, $E$. (c) The number of energy minima increases subexponentially as number of neurons included in the model increases. Error bars are the standard deviation of 10 different subgroups of $N$ neurons. (d) The rank-frequency plot for frequency of visiting each basin matches well between data and model for a typical subgroup of 40 neurons. (e) The number of metabasins, grouped according to the energy barrier, diverges when the energy threshold $\mathrm{\Delta }E$ approaches 1 from above.

Figure 9

The heat capacity is plotted against temperature for models with different number of neurons, $N$. The maximum of the heat capacity approaches the operational temperature of the C. elegans neural system $T_0=1$ from below as $N$ increases. Error bars are the standard error across 10 random subgroups of $N$ neurons.

Figure 10

The topology of the learned maximum entropy model approaches that of the structural connectome, as the number of neurons being modeled, $N$, increases. The two global topological properties being measured are the clustering coefficient $C$ [panel (a)] and characteristic path length $L$ [panel (b)]. Here, the inferred network topology for three different worms is plotted in blue circles. Red crosses are for the randomized network with the same number of neurons, $N$, and number of connections, $N_E$, as the model, where we expect $L_{\text{random}}\sim ln\left(N\right)/ln(2N_E/N)$ and $C_{\text{random}}\sim 2N_E/N^2$. The dark blue line corresponds to the network property of the structural connectome; the dark red line corresponds to randomized network with number of nodes and edges equal to those of the structural connectome [45]. Error bars are generated from the standard deviation across different 10 subgroups of $N$ neurons.

Figure 11

Local perturbation of the neural network leads to global response. (a), (b) For a typical group of $N=50$ neurons, the inferred interaction matrix $J$ is sparse. Here, the neuron indexes $i$ and $j$ are sorted based on $m_i^{\text{flat}}$, as in Fig. 4. (c), (d) When neuron $k$ is clamped to a constant voltage, the Kullback-Leibler divergence (in bits) of the marginal distribution of states for neuron $i$ is distributed throughout the network. (e), (f) When neuron $k$ is ablated, the $D_{\mathrm{KL}}$ is also distributed throughout the network, but is smaller than in response to clamping.

Figure 12

Equilibrium dynamics of the inferred pairwise maximum entropy model fails to capture the neural dynamics of C. elegans. The mean occupancy time of each basin of the energy landscape, $\langle {\tau }_{\alpha }\rangle$, is plotted against the fraction of time the system visits the basin, $P_{\alpha }$. For 10 subgroups of $N=10$ (in dots) and $N=20$ (in asterisks) neurons, the empirical dynamics exhibits a weak power-law relation between $\langle {\tau }_{\alpha }\rangle$ and $P_{\alpha }$. The striped patterns are artifacts due to finite sample size. In contrast, equilibrium dynamics extracted from a Monte Carlo simulation following detailed balance shows an inverse logarithmic relation between $\langle {\tau }_{\alpha }\rangle$ and $P_{\alpha }$, which can be explained by random walks on the energy landscape. Error bars of the data are extracted from random halves of the data. Error bars of the Monte Carlo simulation are calculated using correlation time and standard deviation of the observables.