2.1.

Image Preprocess: Motion Registration, Denoise, and Hemodynamic Correction

Calcium imaging is often accompanied by motion artifacts, even if the animal was head-fixed and anesthetized. In head-fixed experiments, the non-rigid warping of the brain tissue could be caused by heart-beating, breathing, or the shrinking of the tissue from exposure to the laser. Moreover, with freely moving animals, the motion artifact becomes even more severe, because the posture changes induce relative movement of the head and objective lens. In long-term observation experiments, which expand several hours or days, as well as experiments involving different animal subjects, the distortion between frames becomes relatively severer.

We usually assume the stability prior during the calcium extraction step. Thus, an efficient motion-correction algorithm improves the accuracy of neuron extraction. Motion registration is often conducted by transforming each frame toward a reference frame using a mapping function. Linear transformations could be described by a linear matrix. Consider a 2D images as an example (3D volume could be generalized through adding an element to the coordinate vector). The mapping function between the coordinates $x={[x_1,x_2,1]}^T$ and $x^{\prime }={[x_1^{\prime },x_2^{\prime },1]}^T$ of the corresponding sample points could be described as a $3\times 3$ linear transformation matrix $T$, where $x^{\prime }=Tx$. The number of free parameters in $T$ determines the transformation type, which could be divided into rigid transformation (translation and rotation) and non-rigid transformation (similarity, affine, and scaling)⁵⁵ [Fig. 2(a)]. Nonrigid or nonlinear transformations could also be realized by first splitting the image into overlapping patches to perform rigid correction, and merging the patches inversely for an partition-based non-rigid transformation field.⁵⁷ In head-fixed experiments, rigid transformations are usually sufficient for correction of a FOV smaller than 1 mm, and in mesoscale-imaging and freely moving animals, the non-rigid registrations are often demanded. 3D volume non-rigid registration raises the challenge of high computational load. The implementation of the graphic processing unit (GPU)-enhanced algorithms may be considered in large dataset registration missions. In some specific algorithms, the acceleration rate could be up to 100-fold.⁵⁸

Fig. 2

Inter-frame and intra-frame motion artifact correction. (a) Inter-frame motions are global transformations of the whole image which include rigid transformation (translation, rotation, and uniform scaling) and non-rigid transformation (scaling, affine, and projective). (b) Intra-frame motions render as pixel-wise displacement, which is usually caused by the motion during point scanning in multi-photon microscopies. (c) Diagram of intensity-based registration methodologies for inter-frame motion registration. (d) Example estimated trace produced by Lucas and Kanade algorithm, showing the displacement of the specimen in $x$ (black) and $y$ (cyan) directions in each frame and displacement difference between the two frames. Panels (b) and (d) are adapted from Ref. 56. Panel (c) is adapted from Ref. 55.

There are both inter-frame motion artifacts and intra-frame motion artifacts in neural imaging. The former exists widely in almost all kinds of microscopic modalities with the morphological shifting of the brain tissue. And the latter exists mainly in MPM because of the pixel-level temporal incoherence of sampling by point-scanning strategy [Fig. 2(b)]. Inter-frame motion artifacts are relatively more rigid than intra-frame motion artifacts, because there is usually no distortion within a single image, and a reference image is easier to find. A reference frame could be determined by visual inspection or averaging adjacent frames. All other frames could be registered to the reference frame.⁵⁹ The features used for registration could be pixel intensity,⁶⁰ extra structural channel,⁶¹ or exogenous landmark.⁶² A similarity measurements, such as pixel-wise difference, correlation, and information theory-based indexes, can also be used to evaluate the difference between the reference frame and the moving frame. An extra regulation term of deformation was often included, to ensure the sparsity of the deformation field⁶³ and avoid overfitting. Finally, various optimizers, mostly based on gradient descend theory, were designed to search for the best transformation parameter of the mapping function [Fig. 2(c)]. A thorough review of inter-frame motion artifact correction was given by Oliveira and Tavares.⁵⁵ Intra-frame motion artifacts were commonly corrected based on the Lucas and Kanade⁶⁴ algorithm [Fig. 2(d)] or a hidden Markov model,^65,66 where each pixel (or each line) was considered to be sampled from a translated tissue at an independent time-point.^56,67 The Lucas and Kanade algorithm expresses the difference between the registered frame and the template frame as the function of the $x\text{-}y$ trajectory. The optimal estimation of trajectory was derived by setting the first-order Tayler expansion of the error function to zero. And the algorithm iteratively updates the displacement in the $x\text{-}y$ axis, until it converges to an optimal $x\text{-}y$ displacement trace which minimizes the difference to the template. To reduce the computational cost, hierarchical approaches of image registration were proposed, where the imaging video was first decomposed into stable and non-stable sections, and different levels of registrations were assigned to the sectioned images.⁶⁸

Apart from lateral motion, axial drifting eventually drew the researcher’s attention with the evolution of increasing imaging time and resolution demand. In 2D parallel imaging, $z$ axial motion correction could be accomplished by a hardware-implemented strategy of real-time focal plane adjustment,⁶⁹ or by computational approaches of calculating the correlation between the time-sequence $z$ slices and a reference cube.^61,65 In high-speed volumetric imaging, the defocus problem is naturally resolved in a certain axial range because it could capture multi-focus planes at high speed. The $z$ axis displacement could be equivalently resolved like the $x\text{-}y$ displacement by simply adding one more dimension to the registration pipeline. Since the axial displacement could not be avoided completely, volumetric imaging is much more robust to 3D motion artifacts compared with plane-scanning or point-scanning approaches.

A second image processing issue is denoising. Higher-SNR images are always preferred since it enhances the efficiency and fidelity of neuron detection and signal extraction. The simplest solution to acquire a higher-quality image was to use higher laser intensity. But it comes with photobleaching, nonlinear phototoxicity, and heating damage, which interferes with the fundamental neurophysiological phenomenon.⁷⁰ Meanwhile, it is worth noticing that under lower light conditions, photon shot noise becomes comparable to the readout noise of the camera sensor, where the hardware implementation could no longer eradicate noises. Therefore, to facilitate subsequent analysis, data-driven methods become critical.

Image denoising is routinely done by balancing the data fidelity and prior knowledge.⁷¹ Commonly used priors in fluorescence imaging include sparsity prior,⁷² also known as low-rank prior,^73,74 and sensor physics-based noise distribution prior.⁷⁵ Sparsity prior implies that spatial-temporal adjacent blocks in the image would have similar distributions. In other words, high-frequency components are dominated by unwanted noises and should be suppressed selectively. Different methods were used to exploit the sparsity priors. The most general method is to add a regularization term to the image reconstruction target function, which constrains the coefficient density in the transformation domain.⁷⁶ The deconvolution or image reconstruction optimization problem subsequently becomes a multi-target problem. The multiple targets could be decoupled under an alternating direction method of multipliers framework, leading to hybrid iterations of image reconstruction and denoising,^73,77 where the data fidelity term and the sparsity term are optimized iteratively. Within the sparsity sub-step, the parameters of the data fidelity term were fixed and the sparsity prior term was optimized, vice versa in the fidelity sub-step. There are various strategies to meet the sparsity prior. A straightforward approach is setting a threshold (soft or hard) cut-off in transformed domain^72,75,78,79 (Fourier transformation domain, cosine transformation domain, wavelet transformation domain, or a learned over-complete dictionary). The cut-off suppresses the high-frequency components and makes sure that the image frequency spectrum dominantly concentrates on low-frequency bands. Sparsity could also be attained using block-matching approaches. In block-matching-based algorithms, spatial-temporal adjacent blocks within the searching window were vectorized and concatenated into a matrix. The singular value decomposition was performed on the matrix, and then a hard or soft threshold was set to suppress the low-energy components.⁸⁰ It was proved that this kind of strategy could efficiently reduce noise and preserve the details.^72,78

Deep-learning methods have facilitated great advances in bioimage denoising.⁸¹ Artificial neural networks were used to explore the underlying features and recover the noise-degraded images efficiently. There are two kinds of training strategies to train a network, supervised and unsupervised. The supervised strategy requires corresponding noisy and noise-free images to work as training data, and iteratively optimizes the network parameter using gradient descent strategy. In structural imaging in vivo, the training dataset could be acquired by pre-tests, where laser-insensitive specimens are used as subject.⁸² And the pre-trained network could then be used to denoise the laser-sensitive specimen under low light conditions. While in functional imaging, this strategy might be problematic due to the non-repetitiveness of calcium transients, leading to the lack of training data for supervised learning. Therefore, self-supervised training was applied in functional neural image denoising issues. The network used solely noisy data as input and output in training. By utilizing the Noise2Noise framework⁸³ and temporal redundancy prior,⁸⁴ the network successfully produced noise-free images in test sessions. Thus, the self-supervised framework overcomes the obstacle of dataset shortage and is also practical in functional neural image missions.

Resolution maintenance is a vital concern in denoising algorithms designed for microscopy. Overall, the state of art frameworks such as local block-matching strategy and deep-learning network have less resolution loss compared with conventional transformational domain cut-off algorithms. But the performance is usually sample-specific and optical system-specific, and could not be easily summarized. The robustness is particularly a vital problem in deep-learning methods. Actually, there is a risk of resolution loss in every denoising algorithm because imaging denoising is essentially an ill-posed inverse problem. And there is a fundamental trade-off between detail preservation and denoising performance, in other words, the data fidelity and the prior knowledge. The ability of resolution preservation depends mainly on the accuracy of modeling the signal distribution and noise distribution, which requires the expertise on imaging system and sample property from researchers.

The following step of image preprocess is hemodynamic correction, which is typical in wide-field microscopy but is insignificant in two-photon microscopy. The blood flow in vessels has an impact on the background fluctuation of the neural image, which is exhibited as a large-variance and low-frequency background signal component. This component is calcium-irrelevant and should be removed. The most common approach to correct the hemodynamic is adding an extra reference channel with a specific excitation wavelength whose emission is calcium-independent.^16,85 The captured fluorescence intensity was divided by the reference channel, $F_c=\frac{F}{R}$, where $F_c$ is the corrected signal; $F$ is the calcium channel signal and $R$ is the reference channel signal. The corrected normalized signal could be further expressed as $\mathrm{\Delta }{F}_c=\frac{Fc-F_{c0}}{F_{c0}}=\frac{F/R}{F_0/R_0}-1=\frac{\frac{F_0(1+\mathrm{\Delta }F)}{R_0(1+\mathrm{\Delta }R)}}{\frac{F_0}{R_0}}-1=\frac{1+\mathrm{\Delta }F}{1+\mathrm{\Delta }R}-1\approx \mathrm{\Delta }F-\mathrm{\Delta }R$, where subscript 0 stands for averaged signal; $\mathrm{\Delta }F=\frac{F-F_0}{F_0}$, $\mathrm{\Delta }R=\frac{R-R_0}{R_0}$; the last equation is derived through a first order Taylor expansion. Another approach without extra reference channel is modeling the hemodynamic as background signal while performing calcium extraction. This computational approach relies on priors of calcium trace patterns and has no extra hardware cost. Detailed descriptions can be found in the next section.

2.2.

Neuron Segmentation and Calcium Trace Extraction

After the image processing pipeline, the fluorescence images are low-noise and spatially settled and are ready for neuron segmentation and signal extraction. This session describes the methods used to extract fluorescence traces (or equivalent spike trains) and spatial footprints from the raw optical images. The data size should thus be reduced to several MBs which filters out the background fluorescence from neuropils, gliocytes, hemodynamic, and out-of-focus background while keeping all the neural-coding information intact. Compared with the raw optical video which might be up to hundreds of GBs or even TBs, this process could be regarded as data compression or dimensionality reduction. Further neural functional and network inference could entirely rely on the extracted spatial-temporal footprint matrix.

An intuitive way of signal extraction is first segmenting individual neurons as regions of interest (ROI), and then calculating the (weighted) average temporal brightness of these pixels as the temporal trace. This kind of approach relies heavily on the imaging quality and meticulous identification of the shape and size of the neurons, which could be diverse between different bio-sensors and imaging modalities and often requires the expert knowledge of the researcher. Analytical approaches would assume neuron somas to be roughly circular in shape and flicker at a certain frequency periodically. Based on these assumptions, several image segmentation methods are used. First, the pixel-wise maximum deviation from the average brightness^11,86 or the standard deviation (SD) over time is derived to form an active map, from which neuron ROIs could be highlighted. Segmentation methods are subsequently used, ranging from manual approaches^11,65 to automatic algorithms. Manual ROI selection may be laborious but assures high ROI quality, which is suitable for small datasets. There are several software implementations available for manual ROI selections with the aid of automatic initializations, such as ImageJ,⁸⁷ SIMA,⁶⁶ and SamuROI.⁸⁸ Automatic algorithms for the segmentation are mainly based on computer vision theories, such as kernel filtering⁸⁹ and graph-cutting theory.⁶⁶ Deep-learning-based methods have also achieved state-of-art cell-segmentation performance,^90–93 while the conventional shortcomings of lack of training dataset, computational cost, and algorithm robustness to different imaging modalities and SNR levels are gradually increasing. Separation of overlapped neurons is a major challenge currently faced by deep-learning approaches. The anisotropic resolution further aggravated this problem. Different deep-learning methods adapted diverse strategies to tackle this issue. For instance, in U-Net⁹¹ framework, the boundary between cells was artificially inserted into the mask of the training dataset, and the corresponding weight of the ridge was increased to force the network to learn the boundary, so overlapped cells were forced to be split into two non-overlap parts. In STNeuroNet,⁹² overlapped neurons were split using watershed algorithm, and the temporal trace was demixed using a linear regression approach.⁹⁴ This would allow overlapping neurons to be separated spatially and temporally. And in Shallow U-Net Neuron Segmentation (SUNs),⁹⁵ segmentation was done frame-by-frame and followed by a merging procedure, in which case overlapped neurons firing at distinct frames could be separated. In the meantime, recent methods with multi-view imaging and reconstruction⁹⁶ can alleviate the anisotropic resolution, which may also increase the fidelity of segmentations.

The main drawback of the active-map-based segmentation is its high missing rate for inactive neurons and neurons with lower fluorescence indicator expression [Fig. 3(a)]. Another commonly used neural detection approach is matrix factorization. The core framework of these methods is to factorize the spatial-temporal matrix $F\in {\mathbb{R}}^{N\times T}$ into spatial component $S\in {\mathbb{R}}^{N\times K}$ and temporal component $A\in {\mathbb{R}}^{N\times T}$ [Fig. 3(b)]. Each column of $F$ containing $N$ pixels is one of the $T$ vectorized frames of the microscopic video. Ideally, $S$ should contain $K$ neurons’ shapes and locations in every column. And $A$ contains the corresponding temporal traces of these neurons. Additionally, there are also noise $E$ and background signal $B$ in $F$, which could be written as

Fig. 3

Calcium extraction via ROI-based approaches and matrix decomposition. (a) Comparison of ROI analysis and decomposition methods. The ground-truth responses are marked by gray lines. (a1) A simulated image which contains three spatially overlapping signal sources colored in red (Bergmann glial fibers), blue and green (Purkinje cell dendritic trees). The background contains black vessels and bright interneuron somata. (a2) ROI analysis identified spatial filters (left) of cellular components and their temporal traces (right). (a3) Matrix decomposition method reveals the spatial footprint (left) and their temporal traces (right). The signal estimation is more accurate than the ROI based-approach shown in (a2). (b) Schematic illustration of matrix decomposition methods on neural video. The video matrix was factorized into $K$ rank-one matrixes, each of which stands for the location and temporal activity of an individual neuron. (c) The performance comparisons of PCA/ICA, CNMF, and CNMFE on simulated data. Panel (a) is adapted from Ref. 97. Panel (c) is adapted from Ref. 98.

If assuming the spatial location of each neuron does not shift through time, which is facilitated by the previous motion registration process, $SA$ could be further factorized into sum of $K$ neuron’s temporal traces. Thus $SA$ could be represented as the sum of $K$ rank-1 matrixes.

There are infinite solutions to the factorization problem if non-constrains are set. If using only the least square of error as the criterion, the formula leads us to rank-$k$ approximation problem, which is well solved by singular value decomposition, also known as principal component analysis (PCA). However, the PCA method alone might be ill-suited to extract single-neuron signals, since every principal component might contain information from multiple cells.⁹⁷ A further assumption is made about the statistical independence of each component, which leads to independent component analysis (ICA). ICA was proved to outperform PCA in identifying single-neuron cells.⁹⁷ Sparsity prior is another prior apart from independence. Sparse heterarchical matrix factorizaion⁹⁹ with dictionary learning is used to segment neurons and then cluster them into hierarchical functional clusters and reveal the network structure of functional circuits. Moreover, stricter priors of the shape and calcium response pattern could be made¹⁰⁰ which narrows down the solution space furthermore. Non-negative matrix factorization (NMF) method puts non-negative constraints on spatial and temporal matrix,^30,101 which is intuitive considering the non-negative nature of optical images and calcium activity traces. The NMF method performs better on noisy data compared with ICA and was widely adopted by a series of improved algorithms, such as constrained nonnegative matrix factorization (CNMF),¹⁰² CNMF for microendoscope data (CNMF-E),⁹⁸ and CNMF with M-estimator to background¹⁰³ [Fig. 3(c)]. However, CNMF with 3D long-term mesoscale video dataset faces the burden of large computational costs up to thousands of GPU-hours.¹⁰⁴ Several accelerated algorithms are reported, such as seeded iterative demixing (SID)¹⁰⁴ and online deconvolution of calcium image.^105,106 Deep-learning algorithms can also be exploited in the future to further reduce the computational costs with more data priors.³¹

The fluorescence trace normally provides sufficient information as spike trains, since the two modalities are convertible through convolution and deconvolution. However, when precise timing information is wanted, deconvolution may be conducted to improve the data quality.⁵⁴ Under CNMF and OASIS framework, the binary spike train is embedded under the autoregressive model of calcium impulse response and could be achieved through primal optimization. If the segmentation was done using ROI strategies, where spikes were not directly accessible, greedy template fitting was commonly used for spike event detection.^8,107,108

Another common operation on the fluorescence trace is signal normalization, which extracts the mean activity of each neuron and calculates the difference from average score as $\mathrm{\Delta }F/F$. These normalizations could further exclude the influence of incoherent background fluorescence and facilitate the analysis of neural functions.

3.1.

Functional Mapping on Single-Neuron Resolution

An emerging challenge for mesoscale optical recording of neural signals is that task-relevant neurons only take up a small proportion of captured neurons and might be spatially dispersed over the FOV.¹⁰⁹ Most of the neurons over the cortex encode spontaneous and uninstructed movements unrelated to the task.¹¹⁰ Including all these neurons in analyses would cost large computational costs and introduce undesirable noises, which reduce the SNR of the neural encoding space. Thus, a first and foremost step before functional analysis is to discriminate the task-relevant neurons from the massive number of entire neurons.

A forthright intuition is to select neurons that exhibited a more active firing pattern during the task [Fig. 4(a)]. The degree to which task events modulate neuron activities could be evaluated through statistical significance tests by comparing the averaged firing rate between task trials and baseline.¹⁰⁹ Diverse time windows could be used, such as pre-stimulus, post-stimulus, and post-movement onset, depending on the interested task variable. However, the significance test employs nothing on the unique distribution of individual neuron activity prior. Neurons with periodic firing or with extremely active or inactive firing properties may be misjudged by statistical tests. Thus, an alternative option is to perform a shuffle test on the same neuron’s activity trace. Shuffling the start time stamp of behavior segments¹¹ or shuffling the timepoint of calcium events⁸ are both feasible. The same firing rate average analysis is performed on the shuffled data. The significance of task modulation is evaluated through how many percentiles of the shuffled data the primary data could exceed.

Fig. 4

Examples of functional mapping and functional circuit identification. (a) Schematic of task-relevant neuron evaluation. About three example neuron traces are shown with three stimuli as event triggers (dashed lines). In statistical test approaches (top), the fluorescence intensity in the pre-stimulus baseline window (shaded in gray) is compared with post-stimulus time windows (shaded in green), deriving the significance index. In shuffle test approaches (middle), post-stimulus reactions in the shuffled dataset are compared with observed true reaction, where a $P$-value is defined as the percentile of exceeded shuffled data by the observed data. In the kernel regression (bottom) method, the linear model was used to predict the neuron response from the event design matrix. The link between the task and the neuron activity is evaluated via the regression coefficient $\beta$. (b) Example weight maps of event kernel for right visual stimulus and nose movement at pixel-level. (c) Diagram of correlation-based functional clustering at the voxel level. The correlation matrix of the pre-selected supervoxels is manually inspected and a reference trace is selected. Each voxel in the brain is correlated to the reference trace, forming a correlation map. Green, positive correlations; magenta, negative correlations. Panel (b) is adapted from Ref. 110. Panel (c) is adapted from Ref. 60.

Significance tests revealed the overall task modulation on the firing rate during a task period. But within a complicated task, there might be more than one task variable, such as engagement, choice, precision, and body gestures. Statistical tasks based on mean firing rate perform badly on separating these fine variables or revealing the non-linear neural representation. To uncover the relationship between single-neuron activity and animal behaviors, superior data analysis methods should be employed.

Linear regression, also known as general linear modal is an effective method for exploring data correlations. The regression modal could be expressed as $y=X\beta$,¹⁰⁸ where $\beta$ is a concatenated event kernel vector to be regressed. $X\in {\mathbb{R}}^{T\times L}$ is a Toeplitz design matrix concatenated by $K$ events of interest $X_k$, where $X_k\in {\mathbb{R}}^{T\times L_k}$, $\sum _{k=1}^K{L}_k=L$, $X=[X_1,X_2,\dots ,X_K]$. $X_k$ contains diagonal 1s at each event onset and corresponding time lags and 0 elsewhere. $y$ is the normalized neuron trace time course as the prediction output. The link strength between the neuron trace and the behavior variable could be evaluated through the coefficient in $\beta$. A larger absolute value in $\beta$ stands for a higher correlation between neuron and behavior. The spatial-temporal map of $\beta$ could also indicate the information flow between different neuron clusters. To avoid overfitting, different regulation terms are added to rescript the sparsity of the coefficient, which make the model more practical, such as reduced-rank kernel regression,¹⁰⁹ ridge regression¹⁰⁸ or lasso regression¹¹⁰ [Fig. 4(b)].

In perceptual pathways, such as tactile, olfaction, and vision, neural coding mechanisms are often better explored, thus allowing more accurate and complicated modeling of neural activities. For instance, in the tactile sensory pathway, the functional mapping modal assumes that the sensory input goes past static, point-wise nonlinearity, and is convolved with a temporal kernel before being input to a Gaussian noise model.⁸ While in the olfaction pathway, different olfactory is believed to be encoded by population vector instead of single-neuron activities.⁹ And in the visual pathway, the ganglion neurons are well known to be modulated by surrounding neurons, forming the widely recognized receptive field.¹¹¹ These models successfully predict the activity in the local circuit with high accuracy and specificity.

Comparing the modal-free analysis with the model-based analysis, the former is more helpful for rough neuron sorting from mesoscale images, while the latter uncovers the explainable mechanism of neural encoding. Thus, we advocate taking advantage of the high data-throughput property of mesoscale imaging as a neural information database, from which neuron clusters could be distinguished and separated depending on their functions. While the subsequent modeling of neural encoding could be more targeted and prior-based. This may extend the boundary of controversial single-brain-area functional mapping restriction.

3.2.

Brain-Wide Single-Neural Level Functional Network Study

It is widely recognized that different levels of clusters of neurons are recruited for diverse behavior tasks.¹¹² Even a simple perceptual input could arouse distributed neural activities across the brain. Communication across multiple regions is crucial for the brain functioning as a system.^29,109,113 Observing the brain in vivo with behaving animals enables the inference of functional circuits across brain areas, by analyzing statistical dependencies between neurons.

Restricted by the detection capability, conventional network inference either concentrated on a local circuit or utilized low-spatial-resolution brain activities through functional magnetic resonance imaging,^113,114 electroencephalogram, or wide-field microscopy.²⁹ Brain activities acquired through these techniques characterize the area-averaged multi-unit neuron activity containing multi-frequency oscillation components,¹¹⁵ thus are temporally smoother and more stable across trials. Single-neuron activities, however, exhibited larger variance, sparsity, and randomicity through trials.¹¹⁶ The commonly used spectrum methods on field potentials decompose different frequency components and analyze the phase lag between brain areas to infer the information flow through the precedence of brain areas. These methods became meaningless on single-neuron traces because the frequency components lack interpretable physical meanings. Metrics of similarity such as instantaneous phase, phase lag, and phase synchronization index¹¹⁷ are no longer applicable. And the quadratic increase of the computational complexity while calculating pair-wise similarity of the neurons also brings challenges to the network inference. All these barriers bring great difficulties and opportunities to interpret the working mechanisms of neural circuits at the single-neuron level with mesoscale neural imaging.

Despite the inefficiency of phase decomposition-based methods, the covariance or correlation remains the most straightforward model-based metric which describes the undirected similarity between two fluorescence signals. The similarity matrix represents the pairwise correlation among the neural population. With the pairwise correlation which describes the similarity of neurons, subsequent clustering¹¹⁸ methods could be used to assemble the neurons into separated groups.¹¹⁹ Manually selected seed with Pearson correlation matrix has been used to identify the functional circuits across the larval zebrafish brain⁶⁰ [Fig. 4(c)]. There are plentiful other unsupervised clustering methods available, such as $k$-means clustering, hierarchical clustering, barcode analysis, and graph-based analysis. These methods have been used in different experiments. The $k$-means clustering was used to find the representation structure at single-neuron level in the rat orbitofrontal cortex.¹²⁰ The $k$-means method remains computationally efficient even with high dimensional data, but the result of $k$-means algorithms depends heavily on the hyperparameter $k$, thus requiring pretest or the prior knowledge from the researcher. The hierarchical algorithm was used to cluster the retinal ganglion cells and exhibited high interpretability and visualization ability,¹²¹ However, it may encounter severe computational complexity with higher-dimensional data. The barcode analysis required a character barcoding process before categorizing. For instance, the neurons from the whole brain of the zebrafish were barcoded via their response to certain stimuli and clustered into 256 classes.⁸⁶ What may be troublesome with barcoding is that the number of categories increases exponentially with the bit of barcode in binary coding, which restricts its application scope.

Apart from the model-based clustering, functional cell assemblies could also be identified through model-free decomposition analysis. Factorization methods are performed on the 3D tensor stacked by multi-trial 2D neural activity maps, which includes PCA, demixed principal component analysis (demixed PCA),¹²² and non-negative tensor factorization.¹²³ These models factorize the tensor into rank-1 components and extract the temporal and trial dimensional vector as trial-consistent and within-trial latent variables of the neuron representation. Though unsupervised, these vectors may correspond to behavioral observations such as performance accuracy and task engagement.

The temporal correlation-based network mentioned above often induces false-positive connections between indirectly coupled neurons.¹¹² A sparser binary adjacency matrix could be derived by hard-threshold or $k$-nearest neighbors. On the other hand, Granger causality, or transfer entropy analysis, is a temporal-precedency-based metric that reveals directed information flow from one node to another. Its basic assumption is that the history of a precedent neuron should be contributing to predicting the future of the downstream neuron. Granger causality has been used to investigate the unbalanced distribution of information density in the somatosensory cortex of the mouse.¹²⁴ Additionally, there has been open-source algorithm toolboxes¹²⁵ for performing Granger causality analysis. The computational cost of Granger causality algorithm increases quadratically with the number of nodes. Improved Granger causality algorithms¹²⁶ are developed to decrease computational costs and avoid overfitting. Transfer entropy also serves as an equivalent for Granger causality under Gaussian variables¹²⁷ with a much lower computational complexity. With the help of optogenetics,¹²⁸ it is possible to perturbate in vivo neurons and watch the network in dynamic. Modal-based Bayesian network was used to inference the neural spike trains and connectivity with simultaneous optogenetic perturb and electrophysiological recording.¹²⁹