A quantitative model of conserved macroscopic dynamics predicts future motor commands

Abstract

In simple organisms such as Caenorhabditis elegans, whole brain imaging has been performed. Here, we use such recordings to model the nervous system. Our model uses neuronal activity to predict expected time of future motor commands up to 30 s prior to the event. These motor commands control locomotion. Predictions are valid for individuals not used in model construction. The model predicts dwell time statistics, sequences of motor commands and individual neuron activation. To develop this model, we extracted loops spanned by neuronal activity in phase space using novel methodology. The model uses only two variables: the identity of the loop and the phase along it. Current values of these macroscopic variables predict future neuronal activity. Remarkably, our model based on macroscopic variables succeeds despite consistent inter-individual differences in neuronal activation. Thus, our analytical framework reconciles consistent individual differences in neuronal activation with macroscopic dynamics that operate universally across individuals.

Keywords: C. elegans; calcium imaging; computational biology; dynamical systems; individual differences; neuronal dynamics; neuroscience; systems biology; systems modeling.

Figure 1.. Neuronal activity is consistently different in different individuals.

(A) Calcium signals (${\displaystyle \mathrm{\Delta }F/F}$) recorded in one animal for ∼15 min by Kato et al. (2015). Each row represents a single neuron. The top 15 rows (above the red line) correspond to neurons unambiguously identified in all animals (shared neurons). (B) Trace of the AVA neuron colored by behavioral state as defined by Kato et al. (2015). (C) Neuronal activity of representative neurons plotted as a function of behavioral phase (Materials and methods) in a single behavior grouped by animal. Colored solid lines show mean activity for each animal. Shaded regions show 95% confidence intervals. Mean and confidence intervals are computed across multiple cycles of the same behavior in each animal. The dashed black line shows the mean across all cycles of the behavior in all animals. The cross-animal mean for ALA and RIML (Dorsal turn) is a good approximation of activity in each animal individually. In contrast, the cross-animal mean of RIML (Backward locomotion) does not represent any individual. Thus, activity of RIML is consistently different among individuals during backward locomotion. (D) Probabilities that neuronal activity from different individuals was drawn from the same distribution (Materials and methods) computed for each neuron in each locomotor behavior. Activity of most neurons differs consistently among individuals in at least one locomotor behavior. (E) Attempts to decode the onset of backwards locomotion are successful within each animal individually using either all neurons or the 15 shared neurons. This confirms that neuronal activation is stereotyped in each animal. Yet decoding fails across animals. Probability of decoding onset of backwards locomotion in one animal on the basis of neuronal activity averaged across other four animals is indistinguishable from chance. Thus, averaging neuronal activity across individuals disrupts behaviorally relevant information. Box plots show distribution of decoding rates bootstrapped across animals and bi-partitions of the data into training and validation datasets (Materials and methods).

Figure 1—figure supplement 1.. Projection of calcium imaging data onto Temporal PCs.

Time derivative of neuronal activity in each animal is projected onto it’s corresponding first two principal components similar to Kato et al. (2015). The only preprocessing applied to the data was a Gaussian smoothing with sigma of one time step. Colors correspond to the fictive behaviors assigned by Kato et al. (2015). Note that the projection of each individual reveals clear loops in the dynamics of the nervous system. Kato et al. refer to this projection as Temporal PCA. Note that activity from each animal is projected onto its own coordinate system. Also see Figure 1—video 1.

Figure 2.. Simulations of the dynamics faithfully reproduce neuronal activity and behavioral statistics.

(A) Experimentally observed (top) and simulated (bottom) activity of 15 shared neurons plotted as in Figure 1A. (B) Trace of the simulated AVA neuron colored according to the inferred behavioral state. In contrast to Figure 1, color here expresses neuronal activity as z-score computed for each neuron individually. Behavioral states are assigned based on experimentally observed distribution of behaviors for each point in phase space (blue $\mathrm{\to }$ forward locomotion; green $\mathrm{\to }$ backward locomotion) (Materials and methods). Both A and B are plotted on the same time-axes. Because of under-sampling, transitions between forward and backward locomotions are left unassigned (gray). (C) Dwell time distributions for forward locomotion, backward locomotion and backing bouts (blue $\mathrm{\to }$ experimentally observed; orange $\mathrm{\to }$ simulated). Backing bouts were defined as repeated episodes of backing behavior separated by short forward locomotion states lasting at most 30 frames ( 10 s) (e.g. black line in B). Although the manifold is constructed on the basis of transition probabilities between states separated by <1 s., the manifold successfully predicts statistics of behaviors over 100 s.

Figure 2—figure supplement 1.. Neuronal activity trace spectral residues.

Each line gives the difference in fractional power of the spectrum of the observed neuron and it’s simulated partner. The simulated and observed spectra are in good agreement for all frequencies above 0.05 Hz (where finite data artifacts are most prevalent).

Figure 2—figure supplement 2.. Auto correlations and auto mutual information inform expected delay time scales.

The auto correlation (A) and auto mutual information (B) for the experimentally observed neuronal activation in the (Kato et al., 2015) dataset (blue boxes), along with the auto correlation and auto mutual information of their derivatives (red boxes) is plotted for time lags up to 100 s. In the article, we used autocorrelations to estimate an appropriate choice of $\tau$. We selected a delay time based on the minimum (across neurons) auto correlation time constant (10 frames, 4 s). Choice of $\tau$ is not trivial for a multivariate system. However, the results do not strongly depend on this choice so long as the total delay time is ∼ 50 frames (∼ 18 s) Figure 2—figure supplement 3.

Figure 2—figure supplement 3.. Effect of delay embedding parameters on model accuracy.

The delay embedding parameters (number of delays and delay length) are varied to test for robustness. Dots represent a unique pair of parameters. For each pair the method is rerun on all 15 neurons from all five animals in the Kato et al. (2015) dataset and the accuracy of this new model’s simulated motor command dwell time distributions are compared to the observed distributions. Color represents the relative information compared to the original model. Relative information of 1 means that the model constructed using a particular set of delay lengths and number of delays performed as well as the parameter values used in the manuscript (red diamond) (Materials and methods). Note that as long as the total delay length is about 50 time steps the model tends to perform well (about 60% of the original model’s performance).

Figure 2—figure supplement 4.. Model can be constructed on the basis of a single neuron.

The method is applied to data from a single neuron and used to build the dwell times of motor commands as in Figure 3. Relative information of 1 means that the model constructed on a single neuron performed the same as the model based on all 15 neurons (Materials and methods). Note that AVAL, AVER and RIML approximate performance obtained when all 15 neurons are used.

Figure 3.. Position in the phase space determines the expected time of future behavioral transitions.

We explore the region of phase space where most instances of motor commands for backward locomotion are terminated and give rise to commands for forward locomotion. (A) The distribution of times to start of forward locomotion (orange line → mean, orange shading → 95% confidence intervals) is simulated as a function of starting phase. Observed distribution of times to start of forward locomotion in 11 new animals from Nichols et al. (blue line → mean, blue shading → 95% confidence intervals) are calculated for each phase bin in the same region. For each phase bin we also computed the average time since onset of backward locomotion. We then used this average time and the dwell time distribution of backward locomotion to calculate the expected time of termination of backward locomotion (purple line → mean, purple shading → 95% confidence intervals) (Materials and methods). This null hypothesis is grossly inaccurate. (B) Times to start of forward locomotion for both the simulation (orange line → mean, orange shading → 95% confidence intervals, filled circles), and the null model built on dwell times (purple line → mean, purple shading → 95% confidence intervals, empty triangles) are plotted against the time to start of forward locomotion command of the newly observed Nichols et al. animals. Individual points show the average time to start of forward locomotion for each individual phase bin (point color). The simulation (orange) is a good predictor of the newly observed data. In contrast, the null model based on dwell times (purple) consistently underestimates the time until initiation of forward locomotion (p < 0.0001).

Figure 3—figure supplement 1.. Predictions of expected time to switch of behavior are still valid without AVA.

AVAL and AVAR were removed from the five animals in the Kato et al dataset and the Nichols et al dataset used for model validation. The prediction of the expected time to motor command transition was performed as described in the article. The removal of AVA does not dramatically effect the results of the model. Times to start of forward locomotion command for both the simulation (orange line → mean, orange shading → 95% confidence intervals, filled circles), and the null model built on dwell times (purple line → mean, purple shading → 95% confidence intervals, empty triangles) are plotted against the time to start of forward locomotion command of the newly observed Nichols et al. animals. Individual points show the average time to start of forward locomotion for each individual phase bin (point color). Note that the relative information as defined in Figure 2—figure supplement 4 and Figure 2—figure supplement 3 is 0.91 for the model without AVA (Materials and methods). This means that losing information from AVA only reduces the model’s performance by ∼10%.

Figure 4.. Manifolds are invariant across animals.

(A) Manifold constructed for a single animal using all 107 neurons in the recording. Manifold is color coded according to locomotor behavior assigned to each snapshot of neuronal activity by Kato et al. The prevalence of each behavior in each phase bin $\mathrm{\Delta }\mathrm{}\theta$ is encoded in the opacity of the color for each behavior (colors corresponding to each behavior are shown at bottom of the figure). Black arrows represent the direction of flux $d\mathrm{}\theta \mathrm{/}d\mathrm{}t$ around the manifold. DPCA1-3 are the first three principle components of the neuronal activity averaged with respect to $\theta$ and $\alpha$. (B) Manifold constructed using the data from all animals using the 15 shared neurons as common basis. Coloring, projection and arrows are the same as in (A). (C) The data from one animal (lines) is projected onto the manifold constructed from four other animals (surface). The trajectories of the left out animal are well-approximated by the manifold built entirely on the basis of other animals. Behavioral states of the left out animal align well with the distribution of motor commands along the phase of the manifold. This allows us to decode the motor command reliably on the basis of position in manifold space ($\theta$,$\alpha$) across animals. The same manifold and data is shown in two different projections in order to emphasize the similarities between manifold and data.

Figure 4—figure supplement 1.. Projection of calcium imaging data onto DPCAs.

Each animal is projected onto its corresponding first 2 DPCA components after delay embedding using our method from the article. Colors correspond to the fictive behaviors assigned by Kato et al. (2015). Note that the projection of each individual reveals the same two loops that are present in the combined manifold found in the article (Figure 4). The loops are illustrated by the transparent black ovals. Thus, our methodology can be used to uncover loops in neuronal dynamics using ∼100 neurons in each animal. Also see (Figure 4—video 1).

Figure 4—figure supplement 2.. Asymmetric diffusion matrix.

Asymmetric diffusion matrix constructed for a single animal using diffusion kernel as in Equation 13 and normalized such that each row sums to 1 (i.e. a Markov matrix). This matrix was then exponentiated four times to simulate diffusion after ∼1 s. This square matrix has dimensions of the number of time points by the number of time points. Element ($i$, $j$) corresponds to the probability that a system found in state ${\mathbf{𝐃}}_i$ at time 0 will be in state ${\mathbf{𝐃}}_j$ after 1 s where $\mathbf{𝐃}$ corresponds to delay embedded neuronal activity. The matrix has two prominent features: patches and diagonal bands. The former correspond to regions in state space where transition probabilities are approximately symmetric and stochastic processes dominate. The latter corresponds to regions where transitions are asymmetric and deterministic flux plays a dominant role.

Appendix 1—figure 1.. Extraction of manifold from incomplete neuronal activity data.

(A) Observation of a single neuron A) from a two neuron system which produces two distinct modes of oscillation. Activity of $A$ at two time points separated by $\tau$, $A\mathrm{}\mathrm{(}t\mathrm{)}$ and $A\mathrm{}\mathrm{(}t\mathrm{-}\tau \mathrm{)}$ (blue dots), is used as x and y coordinates of a single point in the reconstructed phase plane (blue arrows) obtained by delay embedding. (B) When the activity time series in A is projected onto the plane spanned by $A\mathrm{}\mathrm{(}t\mathrm{)}$ and $A\mathrm{}\mathrm{(}t\mathrm{-}\tau \mathrm{)}$, the reconstructed system trajectory reveals two distinct oscillatory cycles present in the complete system. (C) Schematic of asymmetric diffusion map construction. The arrow indicates the direction of preferred flux. Transition probabilities from the state of the system at time $t$, ${\mathbf{𝐗}}_t$, are computed using a local Gaussian kernel (red $\mathrm{\to }$ high probability; blue $\mathrm{\to }$ low probability). Because this kernel is centered on the next experimentally observed state ${\mathbf{𝐗}}_{t\mathrm{+}\mathrm{1}}$, transition to ${\mathbf{𝐗}}_{t\mathrm{+}\mathrm{1}}$ is most likely. However, transitions to parallel trajectories (red shows high probability of transition) are also possible. (D) Asymmetrical diffusion map $M$ computed for all pairs of states. $M$ is sparse because for each point, transition probability is nonzero just for $k\mathrm{=}\mathrm{12}$ nearest neighbors. Exponentiation of $M$, which shows the evolution of the system after $N$ time steps, highlights the asymmetry. The red dashed line shows the main diagonal illustrating that the parallel bands are off center. Eigen-decomposition of the matrix reveals a linear sum of decay equations. ${\lambda }_i$ is the $i$-th eigenvalue and ${\varphi }_i$ is the $i$-th eigenvector. $c_i$ represent the initial conditions of the system. (E) Manifold constructed for the two neuron system using only neuron A. Manifold construction method (Materials and methods) uses $\varphi$ with the slowest decay time-constant to assign each point in delay embedded space (B) a phase $\theta$ (color of the manifold). $\theta$ identifies position of the system along the cyclic flux. Neuronal activity within a phase bin $\mathrm{\Delta }\mathrm{}\theta$ is averaged such that the center of mass of all points within a bin $\mathrm{\Delta }\mathrm{}\theta$ yields the coordinate of this bin in the delay embedded space. This is used to construct a one-to-one map between a cloud of points in the delay embedded neuronal activity and a phase bin in the manifold space (number of points in each $\mathrm{\Delta }\mathrm{}\theta$ is shown by thickness). Thus, as the system evolves along $\theta$ in the direction given by $d\mathrm{}\theta \mathrm{/}d\mathrm{}t$ (black arrows) it traces a trajectory in the delay embedded neuronal activity. Flux separation, accomplished using cluster analysis, assures that similar phases of different fluxes are not mistakenly grouped together (Materials and methods). (F) Simulations of the manifold in the phase space ($\theta$, $\alpha$) recapitulate neuronal activation. To show this we compare the simulated distribution of time intervals between consecutive minima of A (return times) and those produced by the full system including both neurons (minima shown by black Xs in A).

Appendix 1—figure 1—figure supplement 1.. State space and flux in dynamical systems with noise.

(A) Network diagram of the two neuron toy system explored in the paper. ‘Excitatory’ neuron $A$ synapses (red triangle) onto ‘inhibitory’ neuron $B$ which synapses (blue triangle) back onto neuron $A$. (B) Simulation of the two neuron system with noise. Red and blue traces show activity of neuron $A$ and $B$ respectively. (C) Black arrows represent the flux $\mathrm{(}d\mathrm{}A\mathrm{/}d\mathrm{}t\mathrm{,}d\mathrm{}B\mathrm{/}d\mathrm{}t\mathrm{)}$ for positions in the state space ($A$, $B$). Lines represent trajectories spanned by the state of the system as it evolves in time without noise. Trajectories either converge onto the outer loop (purple) or onto the inner loop (green) over time. These two loops are the cyclic fluxes of the system. (D) Adding noise to the system causes the trajectories to tangle and stochastically switch between the inner and outer cyclic fluxes. Nonetheless, the probability density of points is predominantly concentrated around the two cyclic fluxes and evolves in the same direction as in the system without noise.