(A) Simulated activity of N = 500 neurons organized in A = 5 assemblies. Time sequences of length M = 1000 for all neurons have been generated from the model with λ_μ(0) = 0.08, λ_μ(1) = 0.6 and p_μ = 0.1 for asynchrony, synchrony and activity equal for all assemblies. We used this activity matrix as input to our inference algorithm to recover the five assemblies. (B) (top) log-likelihood over the course of 300 iterations of the Monte Carlo sampler and (bottom) number of transitions per iteration divided by the number of neurons N (different colors correspond to different initializations of neuronal memberships). (C) Assignments of cells to assemblies at four stages of the sampling. Cells initially assigned to O(N) assemblies are progressively grouped until the original assemblies are recovered.

(A) Simulated activity of N = 500 neurons organized in A = 5 assemblies. Time sequences of length M = 1000 for all neurons have been generated from the model with λ_μ(0) = 0.08, λ_μ(1) = 0.6 and p_μ = 0.1 for asynchrony, synchrony and activity equal for all assemblies. We used this activity matrix as input to our inference algorithm to recover the five assemblies. (B) (top) log-likelihood over the course of 300 iterations of the Monte Carlo sampler and (bottom) number of transitions per iteration divided by the number of neurons N (different colors correspond to different initializations of neuronal memberships). (C) Assignments of cells to assemblies at four stages of the sampling. Cells initially assigned to O(N) assemblies are progressively grouped until the original assemblies are recovered.

The ability to identify assemblies in a dataset depends on the model parameters. When synchrony and asynchrony levels are too similar, the recording time might not be sufficient to capture small differences in firing patterns. (A) Raster plot displaying detectable (red) versus non-detectable (blue) regimes of as a function of synchrony and asynchrony parameters of the model (we used 5 randomly generated assemblies for each configuration of synchrony and asynchrony). (B) Raster plot representing the number of neurons reassigned, on average, for every random sample of the Markov chain. The average transition rate can be viewed as an order parameter characterizing the transition between detectable and non-detectable phases.

(A) Performance comparison across levels of asynchrony. Dots correspond to independently generated data sets while solid lines show the average performance for each method over all simulated data. (B) Comparison across number of neurons. (C) Comparison across number of assemblies. (D-F) Standard deviation of the performance across simulated data per parametric condition. Unless specified otherwise, surrogate datasets were generated using 400 neurons and 1000 time frames distributed over 5 assemblies with assembly activity of 5%, synchrony 50% and asynchrony 10%. For k-means we used a number of clusters estimated according to the silhouette method. The number of significant principal components in PCA clustering was obtained by the circular shuffling method whereas for the spectral clustering we used the Newman-Reinart graph-theoretic community detection method [22] (see Materials and methods). The performance was measured according to Eq (27) by comparing each set of assignments with ground truth assignments.

(A) Volumetric 2-photon imaging of both tectal hemispheres showing anterior/posterior (A/P) and right/left (R/L) axis of the zebrafish optic tectum. We monitored activity in the tectum of immobilized zebrafish larvae kept in total darkness. In all experiments we recorded calcium activity for 1 hour throughout 5 planes, 15μm apart with an acquisition frequency of 4.8Hz per volume. (B) Raw images (top) were segmented (bottom) to obtain the temporal dynamics of calcium of thousands of neuron in the tectum (see also S1 Video). (C) Raster plot showing calcium activity of all neurons over 1h recording. (D) Raw fluorescence is described in terms of a Hidden Markov model (HMM) where periods of increased activity are indicated by a binary process ({s}) triggering the onset of calcium transients ({c}). The recorded fluorescence is obtained as the sum of calcium level and a Wiener process to account for the low-frequency baseline modulation ({s}) (see Materials and methods). (E) Example of raw fluorescence and estimation of hidden variables of the HMM using a sequential maximum-likelihood algorithm (see Materials and methods).

(A) Number of neuronal assemblies estimated for each fish. (B) Fraction of neurons assigned to assemblies with 99% confidence. The efficacy of our method in capturing neuronal assemblies is illustrated by separating the activity matrix of assembly neurons from neurons which are not part of a coherent assembly. (C) When sorting the activity matrix by membership, neurons that are part of an assembly generate bands of synchronous activity whereas neurons that are not assigned to an assembly display independent random firing events. (D) Comparison of synchrony, asynchrony, activity and assembly sizes for all experiments. Synchrony and asynchrony correspond to the model parameters λ(1) and λ(0), as the probabilities of a neuron to fire with its assembly (synchronous activation) and independently of it (asynchronous activation). The activity corresponds to the probability of an assembly to be active at any time during the 1 hour recording period. (E) Representative assemblies from single fish with high synchrony and activity display spatial compactness within left or right tectal hemispheres. Labels A/P and R/L indicate anterior/posterior axis and right/left optic tectum respectively. Assembly maps are obtained by projecting in 2D the positions of all member neurons across the five imaging planes.

(A) Raster plot of assembly activity ordered according to subnetwork membership as shown in C. (B) Raster plot of the time correlation matrix (Pearson’s) between all pairs of assemblies averaged across posterior samples and sorted according to subnetwork membership. (C) Subnetwork graphs representing the time correlations among neuronal assemblies (see Materials and methods). Each assembly is a node, shorter edges reflect higher temporal correlations. Note the absence of edges between assemblies located in ipsilateral anterior and posterior tectum (see S3 Video). (D) Scatter plot showing the relationship between time correlation and physical distance between neuronal assembly centers (1pixel ≈ 0.9μm). (E) Fraction of assembly pairs with correlation larger than λ as a function of λ. Ipsilateral pairs of assemblies located at opposite sides in the anterior-posterior axis display a significant reduction in correlation compared to neighboring assembly pairs (anterior-anterior and posterior-posterior). (F) Locations of all assemblies from a single larva (see S2 Video for the assembly distributions in 3D). Each assembly is color-coded according to subnetwork memberships shown in C.

(A) 2D spatial distribution of neurons recorded from the mouse visual cortex from Ref. [19]. The color code represents selected assemblies obtained with our method. (B) Assembly activity matrix. (C) Correlation between the activity of each assembly and (from left to right) the time course of running speed, pupil area and the (binarized) pupil contraction. (D) Correlation and anti-correlation of neuronal assemblies with respect to running speed. (E) Assembly correlated with pupil contraction.

(A) Identification of exponentially decaying activity transients from spike counts using bin size of 0.6s. (B) Binary activity matrix obtained by combining transient events from 242 neurons. (C) Raster plot of the assembly activity. (D) Raster plot of the time correlation matrix between neurons sorted by assembly membership. (E) Time correlation between each assembly activity and the absolute wheel velocity. Boxplots represent the distribution of correlation across posterior samples. (F) Comparison between correlated and anti-correlated assemblies with respect to wheel motion.