System and method for mapping of local brain cell activity and connectivity

Abstract

There is provided a system and method for generating a mapping of local brain cell activity and connectivity. The method including: receiving brain imaging data from image captures of a brain; generating a binary segmentation map using a trained deep learning model, the deep learning model taking the brain imaging data as input; performing registration of the brain imaging data to a brain template; voxelizing the segmentation map and warping the voxelized segmentation map using deformations obtained during registration; determining group-wise heatmaps of neuronal density or cellular density from the voxelized and warped segmentation map; generating a map of clusters of between-group differences in neuronal density or cellular density from the group- wise heatmaps of the voxelized and warped segmentation map; and outputting the map of clusters as a representation of local brain cell activity and connectivity.

Description

SYSTEM AND METHOD FOR MAPPING OF LOCAL BRAIN CELL ACTIVITY AND CONNECTIVITYTECHNICAL FIELD

[0001] The following relates, generally, to brain signal interpretation; and more particularly, to a system and method for generating a mapping of local brain cell activity and connectivity.BACKGROUND

[0002] Mapping neuronal activity and morphology is critical for understanding brain network dynamics underlying behavior and cognition. Advances in microscopy, such as light sheet fluorescence microscopy (LSFM), and tissue clearing techniques such as CLARITY, CUBIC, iDISCO, and SHIELD, have enabled high-fidelity imaging of cell bodies and immediate early genes (lEGs) in intact tissue. These approaches provide useful insights into brain structure and function. However, these imaging and molecular approaches produce exceedingly large (teravoxel scale), complex, multichannel three-dimensional (3D) datasets. These exceedingly large and complex datasets generally introduce substantial challenges for neuroanatomy or neurophysiology analyses.SUMMARY

[0003] In an aspect there is provided a method for generating a mapping of local brain cell activity and connectivity, the method comprising: receiving brain imaging data from image captures of a brain; generating a binary segmentation map using a trained deep learning model, the deep learning model taking the brain imaging data as input; performing registration of the brain imaging data to a brain template; voxelizing the segmentation map and warping the voxelized segmentation map using deformations obtained during registration; determining group-wise heatmaps of neuronal density or cellular density from the voxelized and warped segmentation map; generating a map of clusters of between-group differences in neuronal density or cellular density from the group-wise heatmaps of the voxelized and warped segmentation map; and outputting the map of clusters as a representation of local brain cell activity and connectivity.

[0004] In a particular case of the method, the brain imaging data comprises whole-brain imaging data.

[0005] In another case of the method, the registration comprises performing one or more of: an initial alignment using an initial affine transformation; a rigid registration that allows only translations and rotations to refine an alignment; an affine registration that incorporates one or more of translations, rotations, scaling, and shearing; and a non-rigid, deformable registration using a b-spline symmetric normalization.

[0006] In yet another case of the method, the map of clusters comprises a p-value map.

[0007] In yet another case of the method, the map of clusters is generated using at least one of voxel-wise comparisons, cluster definition and enhancement, and permutation testing.

[0008] In yet another case of the method, voxelizing the segmentation map comprises using a convolution filter and warping the segmentation map comprises using a brain reference atlas using the deformations obtained during registration.

[0009] In yet another case of the method, the heatmaps are determined by subtracting an average of the voxelized and warped segmentation map in each group of the group-wise heatmaps.

[0010] In yet another case of the method, the between-group differences are determined by performing cluster-wise threshold-free cluster enhancement permutation analysis between the groups on the group-wise heatmaps.

[0011] In yet another case of the method, the threshold-free cluster enhancement permutation analysis comprises a group-wise two-way analysis of variance (ANOVA) test.

[0012] In yet another case of the method, the method further comprises outputting auxiliary information about the clusters, the auxiliary information comprising at least one of volume, morphological features, location, and portion of each brain region that is included in each cluster.

[0013] In another aspect, there is provided a system for generating a mapping of local brain cell activity and connectivity, the system comprising one or more processors in communication with a non-transitory storage medium, the non-transitory storage medium comprising instructions for the one or more processors to execute: an input module to receive brain imaging data from image captures of a brain; a segmentation module to generate a binary segmentation map using a trained deep learning model, the deep learning model taking the brain imaging data as input; a registration module to perform registration of the brain imaging data to a brain template, and tovoxelize the segmentation map and warp the voxelized segmentation map using deformations obtained during registration; a statistics module to determine group-wise heatmaps of neuronal density or cellular density from the voxelized and warped segmentation map, and to generate a map of clusters of between-group differences in neuronal density or cellular density from the group-wise heatmaps of the voxelized and warped segmentation map; and an output module to output the map of clusters as a representation of local brain cell activity and connectivity.

[0014] In a particular case of the system, the brain imaging data comprises whole-brain imaging data.

[0015] In another case of the system, the registration module performs the registration by performing one or more of: an initial alignment using an initial affine transformation; a rigid registration that allows only translations and rotations to refine an alignment; an affine registration that incorporates one or more of translations, rotations, scaling, and shearing; and a non-rigid, deformable registration using a b-spline symmetric normalization.

[0016] In yet another case of the system, the map of clusters comprises a p-value map.

[0017] In yet another case of the system, the map of clusters is generated using at least one of voxel-wise comparisons, cluster definition and enhancement, and permutation testing.

[0018] In yet another case of the system, the registration module voxelizes the segmentation map using a convolution filter and warps the segmentation map comprises using a brain reference atlas using the deformations obtained during registration.

[0019] In yet another case of the system, the statistics module determines the heatmaps by subtracting an average of the voxelized and warped segmentation map in each group of the group-wise heatmaps.

[0020] In yet another case of the system, the statistics module determines the between-group differences by performing cluster-wise threshold-free cluster enhancement permutation analysis between the groups on the group-wise heatmaps.

[0021] In yet another case of the system, the threshold-free cluster enhancement permutation analysis comprises a group-wise two-way analysis of variance (ANOVA) test.

[0022] In yet another case of the system, the output module further outputs auxiliary information about the clusters, the auxiliary information comprising at least one of volume, morphological features, location, and portion of each brain region that is included in each cluster.

[0023] These and other aspects are contemplated and described herein. It will be appreciated that the foregoing summary sets out representative aspects of the system and method to assist skilled readers in understanding the following detailed description.DESCRIPTION OF THE DRAWINGS

[0024] A greater understanding of the embodiments will be had with reference to the Figures, in which:

[0025] FIG. 1 shows a conceptual block diagram of a system for generating a mapping of local neuronal, or cellular, activity and connectivity, according to an embodiment;

[0026] FIG. 2 is a flowchart of a method for generating a mapping of local neuronal activity and connectivity, according to an embodiment;

[0027] FIGS. 3A to 3D illustrate an exemplary methodological workflow of the system of FIG. 1 in order to determine an unbiased mapping of local neuronal activity and connectivity;

[0028] FIG. 4 illustrates performance of the system of FIG. 1 in brain-wide segmentation of neuronal cell bodies;

[0029] FIGS. 5A to 5C show region-of-interest-wise evaluation of segmentation of the system of FIG. 1 in segmenting neuronal cell bodies across a whole brain;

[0030] FIGS. 6A to 6E illustrate results of a mouse-related example experiment for mapping local neuronal activity changes underlying food-seeking behavior following cold-stress brain-wide in an atlas-agnostic approach using the system of FIG. 1 ; and

[0031] FIGS. 7A to 7F illustrate results of a locomotion example experiment for brain-wide identification of local neuronal activity changes underlying walking using the system of FIG. 1.DETAILED DESCRIPTION

[0032] For simplicity and clarity of illustration, where considered appropriate, reference numerals may be repeated among the Figures to indicate corresponding or analogous elements. In addition, numerous specific details are set forth in order to provide a thorough understanding of the embodiments described herein. However, it will be understood by those of ordinary skill in the art that the embodiments described herein may be practised without these specific details. In other instances, well-known methods, procedures and components have not been described in detail so as not to obscure the embodiments described herein. Also, the description is not to be considered as limiting the scope of the embodiments described herein.

[0033] Various terms used throughout the present description may be read and understood as follows, unless the context indicates otherwise: “or” as used throughout is inclusive, as though written “and / or”; singular articles and pronouns as used throughout include their plural forms, and vice versa; similarly, gendered pronouns include their counterpart pronouns so that pronouns should not be understood as limiting anything described herein to use, implementation, performance, etc. by a single gender. Further definitions for terms may be set out herein; these may apply to prior and subsequent instances of those terms, as will be understood from a reading of the present description.

[0034] Any module, unit, component, server, computer, terminal or device exemplified herein that executes instructions may include or otherwise have access to computer readable media such as storage media, computer storage media, or data storage devices (removable and / or nonremovable). Computer storage media may include volatile and non-volatile, removable and nonremovable media implemented in any method or technology for storage of information, such as computer readable instructions, data structures, program modules, or other data. Examples of computer storage media include RAM, ROM, EEPROM, flash memory or other memory technology. Any such computer storage media may be part of the device or accessible or connectable thereto. Further, unless the context clearly indicates otherwise, any processor or controller set out herein may be implemented as a singular processor or as a plurality of processors. The plurality of processors may be arrayed or distributed, and any processing function referred to herein may be carried out by one or by a plurality of processors, even though a single processor may be exemplified. Any method, application or module herein described may be implemented using computer readable / executable instructions that may be stored or otherwise held by such computer readable media and executed by the one or more processors.

[0035] To enable automated brain-wide activity mapping with immediate early genes (lEGs) in large microscopy datasets, some pipelines rely on registration to standardized brain atlases or common reference spaces for statistical analyses, via a region of interest (ROI)-based approach. These types of analyses require a priori knowledge or data-specific expertise to choose ROIs and compare cell counts across groups. Notwithstanding, emerging single-cell and spatial-omics data and transcriptomic analyses indicate that there is far greater diversity within conventionally defined atlas regions, with neuronal subpopulations with unique cytoarchitecture, connectivity, and function. Moreover, definitions of atlas regions are commonly structure-centric and limited in their delineation of neuronal subtype-specific and local areas. Thus, relying on a traditional regionspecific grouping of voxels can fail to detect subtle focal contrasts, neuronal subpopulationspecific effects, or brain-wide changes, in an unbiased fashion. Further, aggregating results on a regional basis can obscure the heterogeneity of changes within brain regions, yet many pathologies are thought to exert salient laminar or columnar effects. In addition to atlas-based regional analyses, other pipelines for light sheet fluorescence microscopy (LSFM) data enable voxel-wise statistical analyses (e.g. using an independent two-sample Student’s t-test) to assess neuronal changes in different brain regions. However, in-vivo neuroimaging studies have demonstrated that voxel-wise statistical tests with leniently corrected p-values can result in an inflated rate of false positives. This issue is greatly amplified with tera-voxel scale whole-brain LSFM images containing trillions of voxels. Conversely, conservative corrections, such as Bonferroni or Benjamini and Houchberg (false discovery rate) drastically reduce the power of these voxel-wise approaches to detect salient changes in these tera-voxel datasets.

[0036] Commonly employed pipelines for cellular segmentation using fluorescent microscopy images rely on traditional image processing techniques in which parameter tuning or expert intervention is often required to extract meaningful features for segmentation. This limits their ability to produce robust results on unseen datasets, with varying signal and noise distributions across different microscopes, imaging protocols, artifacts, samples, fluorescent reporters, and brain regions. Deep learning (DL), especially Convolutional Neural Networks (CNNs) and vision transformers (ViT), can automatically learn effective representations of data with multiple levels of abstraction, resulting in accurate and robust segmentation of imaging data. While pipelines have been attempted that leverage DL for mapping cells in microscopy data, they are substantially confounded by reliance on 2D-based DL models, which impede the unbiased assessment of 3D volumetric changes. Additionally, such pipelines are substantially limited by training and testing on restricted datasets consisting of a small sample from specific brain regions and cell types;limiting their ability to segment a large variety of sizes, shapes, and densities in different brain regions and distributions. Further, such pipelines are substantially limited by combining conventional image processing techniques with a DL-based classifier, without developing an end- to-end segmentation model; limiting their functionality to detection. Such DL-based pipelines are specifically tailored for 3D mapping of neuronal activity in whole-brain LSFM data, thus hindering their effectiveness in handling the variations across datasets. Furthermore, current segmentation pipelines do not provide uncertainty estimates of model predictions, which are invaluable for evaluating the reliability of segmentation models, and for guiding the improvement of segmentation results.

[0037] The present embodiments advantageously overcome, at least, the above limitations in the art and enable robust brain-wide mapping of neural subpopulation-specific effects in LSFM data. The present embodiments can provide an end-to-end automated pipeline to enable an unbiased and generalizable mapping of the 3D activity or morphometrical changes of neurons in tera-voxel scale LSFM data. Unlike approaches, the present embodiments, via integration with multimodal image registration and connectivity analysis (informally referred to as “MIRACL”) provides a quantitative mapping of neuronal subtype-specific effects to brain regions and clusters that are activated in response to experimental manipulations in an atlas-agnostic manner; i.e. independent of pre-defined atlas regions. The present embodiments advantageously provide threshold-free cluster-wise permutation analysis, optimized for LSFM whole-brain data, enabling cluster-based statistical analysis, which improves the sensitivity and specificity of voxel-wise analysis in LSFM.

[0038] Embodiments of the present disclosure (informally referred to as “ACE”) combine deeplearning (DL), registration techniques, and statistical analysis to map neuronal activity in wholebrain tera-voxel light sheet fluorescence microscopy (LSFM) data in order to assess cluster-wise and neuronal subtype-specific changes. Considering the significance and inherent technical complexities entailed in segmenting neuronal cell bodies across the brain and samples from different microscopes and imaging protocols, deep-learning architectures were used by training them on large datasets of whole-brain LSFM data from different centers. Advantageously, vision transformers (ViT) DL models can be used, in embodiments of the present disclosure, as a segmentation core, providing a quantitative estimate of model confidence (or uncertainty), and cluster / voxel-wise statistical analyses can be performed via an advanced non-parametric permutation-based algorithms.

[0039] Embodiments of the present disclosure can use optimized ViT as a backbone architecture for a deep learning model, which can use a self-attention-based transformer operating on vectors from image patches to model long-range relationships and dependencies in the data, and improve learning of global information. In some cases, the output can be obtained by ensembling a plurality of models’ predictions (e.g., 50 models) using, for example, the Monte-Carlo (MC) dropout technique to estimate model uncertainty (via computing the variance across models) and enhance robustness through averaging predictions generated from stochastically different models.

[0040] Due to its modular architecture, embodiments of the present disclosure are compatible with a wide range of deep learning-based segmentation models. For instance, a trained ACE model designed for segmenting neuronal activity can be substituted with an alternative deep learning model trained specifically for segmenting other cell types in light sheet fluorescence microscopy (LSFM) datasets, without requiring changes to the overall pipeline.

[0041] Embodiments of the present disclosure may, in some cases, employ a student-teacher learning strategy to enhance segmentation module performance. In this approach, pretrained models are utilized to generate silver-standard ground truth labels across various datasets, including application-specific datasets. These generated labels can then be used to further retrain the model, thereby improving its generalizability and accuracy.

[0042] Embodiments of the present disclosure may, in further cases, incorporate a selfsupervised learning approach to enhance the generalizability of the pretrained models. In this approach, an auxiliary task such as image masking and reconstruction is applied to LSFM data that lacks ground truth annotations, enabling the models to learn meaningful representations during pretraining. These pretrained models can subsequently be fine-tuned on datasets with available ground truth to improve task-specific performance.

[0043] In an example experiment, to train and evaluate the DL models and ensure a robust segmentation across brain regions and among heterogeneous samples, LSFM data was acquired from 18 Tg TRAP2-Ai9 mice with cFos promoter. The experiments divided each training dataset into smaller 3D image cubes or patches (963voxels corresponding to 0.353mm3). To generate ground truth (GT) data for each patch, the experiments used a semi-automatic pipeline. Thus, the example experiments used randomly selected 15,200 unique input patches for training, not accounting for data augmentation (N=30,400 patches with augmentation). The example experiments tested an optimized 3D U-Net architecture with residual blocks, and dropout layersas a baseline model to compare. The example experiments illustrated the robustness of the segmentation of the present embodiments on the test set, as well as an unseen dataset, consisting of three LSFM whole-brain datasets acquired from a mouse model of spinal cord injury.

[0044] Turning to FIG. 1 , shown is a conceptual diagram of a system 100 for for generating a mapping of local brain cell activity and connectivity, according to an embodiment. The system 100 can be colloquially referred to as ‘ACE’. As shown, the system 100 has a number of physical and logical components, including a processing unit (“Pll”) 160, data storage 164, an interface module 168, a network module 176, and a local bus 184 enabling Pll 160 to communicate with the other components. Pll 160 can include one or more central processing units, one or more graphical processing unit, microprocessors, dedicated hardware, or other integrated processing circuits. The data storage 164 provides storage of data to the Pll 160 and can include volatile and nonvolatile storage for storing data required by the Pll 160, including computer-executable instructions for implementing the methods described herein, as well as any derivative or related data. The interface module 168 enables input to be provided; for example, directly via a user input device, or indirectly, for example, via an external device. The interface module 168 also enables output to be provided; for example, directly via a user display, or indirectly, for example, sent over the network module 176. The network module 176 permits communication with other systems or computing devices; for example, over a local area network or over the Internet.

[0045] In other embodiments, any operating system, programs, or instructions can be executed in hardware, specialized microprocessors, logic arrays, or the like. While FIG. 1 illustrates a system implemented on a single computing device, it is understood that the processing, or any of the functions undertaken by the system 100, can be distributed over multiple computing devices; for example, in a cloud or distributed computing environment.

[0046] In an embodiment, the Pll 160 can be configured to execute a number of conceptual modules 101 ; for example, an input module 102, a segmentation module 104, a registration module 106, a statistics module 108, and an output module 110. In further cases, functions of the above modules can be combined or executed on other modules. In some cases, functions of the above modules can be executed on remote computing devices, such as centralized servers and cloud computing resources communicating over the network module 176.

[0047] FIG. 2 illustrates a flowchart of a method 200 for generating a mapping of local brain cell activity and connectivity, according to an embodiment. At block 202, the input module 102receives brain imaging data, such as light sheet fluorescence microscopy (LSFM) data, from the data storage 164, the interface module 168, and / or the network module 176. In most cases, the brain imaging data, such as the LSFM data, are images of raw whole-brains that are immunolabeled and cleared; however, in other cases portions of the brain can be received.

[0048] At block 203, the input module 102 performs preprocessing on the brain imaging data (for example, on inputted large-scale whole-brain datasets) to improve computational efficiency. In an example of preprocessing, a three-dimensional (3D) brain mask can be automatically generated using the Otsu thresholding approach to identify and exclude empty regions from further processing. In such cases, these non-informative (or background) patches can be bypassed during deep learning deployment, thereby accelerating computation of the present method. The locations of the excluded patches can be tracked to ensure proper reconstruction of the full 3D segmentation maps.

[0049] At block 204, the segmentation module 104 receives the brain imaging data and generates whole-brain binary segmentation maps using deep-learning (DL) models, such as ViT-based and / or convolutional neural network (CNN)-based machine learning models.

[0050] Using the DL models, the segmentation module 104 can also generate voxel-wise uncertainty maps to estimate the confidence of the DL models. The generated voxel-wise uncertainty maps enable users to determine whether to proceed with the analysis with the deep learning models or to fine-tune these models using the fine-tuning provided herein to accommodate datasets with different image characteristics with respect to their training datasets. Suitable fine-tuning scripts can be integrated to facilitate such approaches. Once fine-tuned, the updated models can be replaced with the default models, and segmentation can be executed accordingly.

[0051] At block 206, the registration module 106 receives the brain imaging data (in some cases, autofluorescence LSFM data) and registers the data to a brain template, for example, the Allen Mouse Brain Reference Atlas (ARA). The registration can include registering the brain imaging data to a common coordinate system via linear and / or non-linear transformations.

[0052] Registration by the registration module 106 brings the whole-brain tissue-cleared microscopy images, segmentation maps, and an atlas into spatial correspondence. In a nonlimiting example, registration can begin by having the microscopy images undergo preprocessing, such as denoising to reduce noise and intensity correction to standardize the intensity values. Inthe example, the registration module 106 can perform a first alignment by utilizing a rough alignment with calculation of an initial affine transformation. This involves estimating translations, rotations, and scaling parameters to approximately position the images in the same spatial space and orientation. The registration module 106 can follow this initial affine transformation by a multistage registration, with the first being a rigid registration with, for example, 6 degrees of freedom; allowing only translations and rotations to refine the alignment. In this example, the registration module 106 can then perform an affine registration with, for example, 12 degrees of freedom; incorporating translations, rotations, scaling, and shearing for further refinement. In this example, the registration module 106 can then perform a non-rigid, deformable registration using a b-spline symmetric normalization to account for complex deformations to achieve an accurate alignment. During the rigid and affine stages, a normalized mutual information similarity metric can be used to maximize the statistical dependence between the images and the reference atlas, while the cross-correlation metric can be employed during the deformable stage to ensure image similarity. The resulting deformations enable bidirectional warping of images, aligning them between the microscopy native space and the atlas space. Registration ensures spatial correspondence, facilitating precise comparative analyses.

[0053] In some cases, the registration module 106 can generate checkpoint visualizations that allow users to assess registration performance and adjust parameters as needed. In such cases, the downsampled input data is overlaid onto the brain template’s annotations, and an animated visualization is generated and provided for user review. The registration module 106 supports re- execution with modified registration parameters independently, without requiring reprocessing of the segmentation step.

[0054] At block 208, the registration module 106 voxelizes the segmentation maps using a convolution filter and warps them to the ARA (10 pm) using deformations obtained from registration. Warping of the segmentation maps can be to a common coordinate framework to enable measurement of neuronal density and morphometric changes within different significant clusters or regions across animals. Warping the maps using deformations obtained from registration generally involves applying the determined transformations to adjust and align the images or maps to match a reference space. These deformations are determined from the registration, which identifies how different voxels of the images need to move to achieve optimal alignment. By warping the maps, the system 100 ensures that anatomical structures are correctly positioned relative to each other.

[0055] At block 210, the statistics module 108 receives the voxelized and warped segmentation maps and determines group-wise heatmaps of neuronal density. The heatmaps are obtained by subtracting an average of the voxelized and warped segmentation maps in each group to identify neural (or cellular) activity hotspots.

[0056] At block 212, the statistics module 108 determines between-group differences in neuronal density by performing region-of-interest (ROI)-wise analysis between the groups on the determined heatmaps. To identify significant localized group-wise differences in neuronal activity in an atlas-agnostic manner, the statistics module 108 can perform a cluster-wise threshold-free cluster enhancement permutation analysis (for example, using a group-wise ANOVA). The analysis generates a p-value map with clusters that represent significant differences between groups. The clusters represent a mapping of local neuronal activity and connectivity. While a group-wise two-way analysis of variance (ANOVA) test is generally preferred, any suitable statistical test / analysis can be used.

[0057] At block 214, the output module 110 outputs the determined clusters to the data storage 164, the interface module 168, and / or the network module 176, the clusters representing the mapping of local neuronal (or cellular) activity and connectivity. In some cases, the output module 110 can output a table that summarizes these clusters, including their volumes and the portion of each brain region that is included in each cluster. In some cases, the output module 110 can output the determined clusters in the form of a visualization; where significant correlation coefficients are plotted alongside a heatmap that illustrates a Euclidean distance between each pair of cluster centroids. In some cases, the output module 110 can output the determined clusters along with auxiliary information about such clusters, such as volume, morphological features, location, and portion of each brain region that is included in each cluster.

[0058] The determined clusters, which can be represented by the p-value map, indicate regions in the brain where statistically significant differences in activity or connectivity between groups have been detected. These clusters, which can be obtained, for example, via voxel-wise comparisons, cluster definition and enhancement, and permutation testing, ensure that the findings are robust and not due to random chance. The present inventors determined that these clusters are indicative of significant group differences in the expression of immediate early genes (or c-fos labeling), which represent neuronal activity after a certain behavior or task. Thus, the statistical strength and spatial extent of these clusters suggest areas of local neuronal activity.

[0059] In some cases, to further assess connectivity, the statistics module 108 can also perform cluster-wise connectivity analysis. The statistics module 108 ranks significant clusters and determines mean intensity values within each cluster for each subject. Pearson correlation tests are can be performed between clusters to determine the relationship and connectivity strength. The output module 110 can output significant correlations identified by the statistics module 108.

[0060] In some cases, to further validate the clusters, the statistics module 108 can also perform native space cluster validation. The statistics module 108 binarizes and differentiates the significant clusters through connected component analysis. These clusters can then be warped back to the native space of each subject using registration transformations. Neuronal counts within each cluster can be determined and comparisons between groups can be determined using statistical tests (e.g., Mann-Whitney II test). The output module 110 can output these counts and / or comparisons.

[0061] While the method 200 generally describes performing the operations of whole brain data, it is understood that the operations can likewise be performed on only portions of the brain instead of the whole brain.

[0062] While the present embodiments generally describe the use of LSFM data as the brain imaging data, it is understood that the system 100 can use any suitable brain imaging modality; for example, multi-photon microscopy data and 3D histology stacks. In this way, the system 100 can be extended to other markers of brain activity or cellular activity in general. Advantageously, the fine-tuning described herein enables adaptation to new datasets with differing characteristics, for example, to identify the targets of a covalent drug in the brain.

[0063] Instead of focusing on individual voxels, the system 100 instead considers a null hypothesis regarding the sizes of clusters in the data by incorporating the correlation structure of the data. Advantageously, clusters are determined in a threshold-free approach that enables the detection of significant changes in regions and sub-regions with higher sensitivity. In this way, the volume, strength (effect size), and the ARA regions spanned by each cluster can be determined. Moreover, integrating cluster-wise results and correlation analysis, the associations and interactions among neuronal ensembles differentially activated in response to experimental manipulations can be determined; potentially revealing both within-region and long-range connectivity or functional coherence.

[0064] FIGS. 3A to 3D illustrate an exemplary methodological workflow of the system 100 in order to determine an unbiased mapping of local neuronal activity and connectivity.

[0065] The present inventors conducted example experiments to verify the substantial advantages of the system 100. The example experiments first evaluated the performance of an ensemble model with the system 100 (using MC dropout and N = 50 models) on a test dataset (consisting of 12,160 unique patches with 963voxels 0.353mm from 5 animals). The high fidelity of the system’s predictions can be visualized FIGS. 3A to 5C on a representative dataset, along with the corresponding GT. To quantitatively assess the performance of the segmentation models, a series of overlap and surface-based metrics were used, including recall, precision, Dice similarity coefficient (DSC), and 95% Hausdorff distance (HD95).

[0066] To provide a fair comparison, llastik (random forest classifier) model was trained using all training subjects. On the test set, the system 100 outperformed llastik across all experiments. The system 100 achieved an average improvement in DSC of 0.17 compared to the optimized llastik model (p < 0.0001). The system 100 showed superior performance in detecting the boundary and shape of the neurons and discriminating neurons that were close to each other, resulting in a lower HD95 (ACE: 4.76±3.47; llastik: 9.60±6.78; p < 0.0001). Furthermore, the system 100 exhibited increased robustness in segmenting neurons with low signal intensity or slight out-of- focus blurriness, while llastik struggled to segment them effectively. The system 100 demonstrated robust segmentation accuracy across different ROIs and consistently superior segmentation performance (on all evaluation metrics) compared to llastik across the whole brain.

[0067] To further validate the segmentation model, the example experiments compared its performance against a detection algorithm called Cellfinder. To this end, the segmentation maps were transformed into detection maps by finding the center of mass of each segmented neuron. The system 100 exhibited a superior detection performance, resulting in an F1 -score of 0.75±0.08 vs. 0.55±0.15 for Cellfinder (p < 0.0001). To improve the accuracy of the Cellfinder pipeline on the experiment dataset, the example experiments retrained (fine-tuned) the Cellfinder’s DL model on the training data and repeated the evaluation experiment. While the retraining elicited an increase of 8.64% in precision compared to the Cellfinder pre-trained model, it did not improve the overall performance, yielding an overall decrease in the F1 -score (to 0.28±0.12).

[0068] Employing an ensemble of 50 ViTs via Monte-Carlo (MC) dropout improved the baseline model performance, resulting in an average of 2.1% improvement in precision. The system’s 100ensemble models showed high confidence (low uncertainty) in correctly segmented regions, whereas mis-segmented regions demonstrated low confidence (high uncertainty). Areas with high uncertainty were typically observed around the boundary of neuronal somas and sparsely miss- labeled processes, such as discontinuous axons. Voxel-wise uncertainty maps provided by the system 100 provide an estimate of segmentation confidence and can be used to remove potentially false positive voxels from prediction maps.

[0069] To validate the generalizability of the DL segmentation models, the example experiments obtained an unseen (out-of-distribution) dataset with different image resolution, scanning parameters, microscope, signal-to-noise ratio, and morphological appearance of neurons from those in the training dataset. From this unseen dataset, 1 ,820 patches of 963voxels each (0.27x0.27x0.48 mm3) were randomly selected and ground truth (GT) data for each image patch was generated. The ensemble model was trained on the training data set and deployed on the unseen dataset without any fine-tuning or post-processing. The average DSC achieved by the system 100 was 0.73±0.02 vs. 0.45±0.12 for llastik (p < 0.0001). Ilastik segmentations resulted in substantial number of false negatives (Recall: 0.35±0.15 vs. 0.78±0.09 for ACE, p < 0.0001). Furthermore, the system 100 outperformed Cellfinder on this unseen dataset (with both the pretrained and fine-tuned models). To ensure a fair comparison, Cellfinder runs were conducted using various parameters for its detection step. From these runs, the experiments selected the best model based on a visual comparison. This approach accounted for the differences in neuronal size distribution and choose the most suitable Cellfinder model for these new datasets. Notwithstanding, the system 100 outperformed the tuned Cellfinder model on all evaluation metrics (p < 0.0001).

[0070] To further increase robustness on unseen data, an additional layer of ensembling was used by combining UNet Transformer (LINETR) and ll-Net ensemble models, generating an ‘ensemble of ensembles’. Segmentation thereby takes advantage of a CNN-based architecture that extracts local features within the image patch and a ViT architecture to learn long-range dependencies across the patch. The ensemble of ensembles strategy increased the generalizability of segmentation when dealing with the unseen dataset, improving the DSC by an average of 2.1 %, precision by 5.2%, and decreasing the average HD95 by 15.4%. In detection mode, the ensemble of ensembles increased the F1 -score by an average of 5.1%, and precision by 11.1%. The consistently superior results support the use of the system’s 100 segmentation for robust and generalizable mapping of neuronal cell bodies across disparate subjects and experiments.

[0071] FIG. 4 illustrates performance of the system 100 in the example experiments in brain-wide segmentation of neuronal cell bodies. In section a, maximum intensity projection rendering of whole-brain c-Fos expression with an enlarged view of a cortical patch is shown. The segmentation maps predicted by the ensemble for the enlarged subregion are shown and compared with ground truth. In section b, the raw image, ground truth, and segmentation maps for two example image patches, along with voxel-wise uncertainty maps are shown. The high uncertainty regions are localized around the boundary of sparsely mislabeled processes such as axons (first column) and neuronal somas (second column). The arrows point to the missegmented regions. In section c, shown is qualitative evaluation of segmentation accuracy of the system 100 vs. Ilastik in terms of detecting neurons with low signal intensity or slight blurriness (first row) and the boundary and shape of neurons (second row). The arrows point to the boundary of two neurons that are close to each other. Section d and e show quantitative evaluation of the segmentation accuracy of the system 100 versus Ilastik and detection accuracy of the system 100 versus Cellfinder in terms of average Dice coefficient, precision, recall, 95% Hausdorff distance, and F1 score on test datasets (N: 12,160 unique patches with 9630.35 mm3) and unseen datasets (N: 1 ,824 unique patches of 9630.27 x 0.27 x 0.48 mm3).

[0072] FIGS. 5A to 5C show ROI-wise evaluation of the segmentation module 104 in segmenting neuronal cell bodies across the whole brain. FIG. 5A shows qualitative evaluation of DL segmentation module in different cortical regions in an example subject from the test dataset. Each panel from left to right demonstrates a 3D maximum intensity projection of a raw (input) image patch, ground truth, model output, and an overlaid version of all three. FIG. 5B shows, on the test set, each LSFM dataset registered to the ARA. ARA labels were then warped to each subject’s native space. The warped labels were used to determine the location of each image patch in the brain. Average Dice coefficient (b) and 95% Hausdorff distance (c) obtained between the outputs and ground truth per ARA label. FIG. 5C shows these locations compared against Ilastik.

[0073] Understanding the neural mechanisms governing cold-induced food-seeking is crucial for unraveling the intricate interplay between environmental stimuli, energy expenditure, and feeding behavior in mammals. The example experiments employed whole-brain c-Fos screening (via SHIELD tissue clearing and LSFM) of mice after a prolonged (6-hour) exposure to either 4°C and 30°C temperatures (n = 4 / group) to map the unique neuronal ensembles that control cold- induced food seeking. To validate the system’s 100 ability to precisely detect selective activationof specific neuronal populations at the sub-regional level in an unbiased manner, the experiments processed and analyzed this whole-brain c-Fos LSFM dataset.

[0074] In the example experiments on the mice, the DL model was used to generate whole-brain segmentation maps from the raw LSFM imaging data. To maintain the fidelity of data and prevent information loss during down-sampling and warping into a common reference with lower resolution, a convolution-based voxelization was used on the whole-brain segmentation maps. Results obtained from automated segmentation of the brain regions were confirmed by inspection of both voxelized segmentation maps and raw data. Following voxelization, segmentation maps were aligned with the Allen reference atlas at a 10 pm resolution using deformation fields obtained. To identify neural activity hotspots and determine imaging signatures reflecting focal cellular changes after cold exposure, group-wise heatmaps of neuronal density were generated by subtracting the average of the voxelized-warped segmentation maps in each group. These heatmaps revealed that cold stress elicited increased neuronal activity within the hypothalamus, consistent with its role in thermoregulation, and attenuated neuronal activity across the cortex, likely due to the reduced physical activity during cold exposure.

[0075] In the example experiments on the mice, to assess between-group differences in neuronal density, a ROI-wise analysis was performed between the two groups within each ARA region. This whole-brain analysis highlighted several brain regions that exhibited significant differential activation (p < 0.05), notably the Paraventricular Hypothalamic Nucleus (PVH) and the Nucleus of Reunions (RE). Some areas showed diverging group trends but did not reach statistical significance, such as the Cuneiform (CUN), a region involved in locomotion, and Primary auditory cortex (AUDp). A visual inspection of the group-difference heatmap revealed that the ROI-based analysis on pre-defined atlas regions failed to detect several highly localized areas with substantial increases in c-Fos activation in the cold group: most notably, within the large RE region. To test whether the system 100 could map significant localized group-wise differences in neuronal activity in an atlas-agnostic manner, a cluster-wise TFCE permutation test was performed using a group-wise ANOVA. The system 100 extracted clusters of activation across the whole brain and summarized and ranked each cluster by its strength of activation (sorted by effect size and p-value), total volume, and its percent volume within each overlapping brain region. Expected areas of c-Fos activation were highlighted and numerous sub-regional changes in the cold-induced dataset were identified, for instance, clusters that were only located within the dorsal sub-region of the paraventricular nucleus of the thalamus; which were missed in the ROI-based analysis using the (whole) PVT region (ARA label). Interestingly, in contrast to the ROI-basedanalysis, the cluster-wise test revealed a localized increase in c-Fos activation in the ventral midline of the RE and above the third ventricle, corresponding to the Xi, and confirming its role in the neural orchestration of cold-induced food-seeking behavior and energy homeostasis. Connectivity analysis was also performed between significant clusters, assessing the correlation of their activations (using Pearson correlation analysis with permutation and bootstrapping). Using the connectivity analysis, significant associations (p < 0.01) were determined between clusters in the midline group of the dorsal thalamus (such as PVT and Xi) and the nucleus accumbens (ACB), ie. putative regional connectivity, which was validated by anterograde viral tracing and the Allen connectivity atlas.

[0076] These results underscore the versatility of the system 100, showcasing its ability to accurately map localized changes in cell activation throughout the entire brain, without reliance on pre-defined atlases; and to generate connectivity maps by correlating significant clusters of activations. This functionality can be further used to provide insights on dynamic networks governing specific behaviors.

[0077] FIGS. 6A to 6E illustrate the mouse-related example experiments for mapping local neuronal activity changes underlying food-seeking behavior following cold-stress brain-wide in an atlas-agnostic approach. FIG. 6A shows an overview of the experiment designed to analyze c- Fos+ cell distribution in whole-brain LSFM data after cold stress (n = 4 / group). FIG. 6B shows that using registration algorithms, whole brain LSFM data were registered to ARA. The left and right upper panels show autofluorescence LSFM data overlaid on ARA labels (boundaries) for two example subjects in each group, showcasing the registration accuracy. FIG. 6B also shows whole-brain segmentation maps of c-Fos+ cells obtained from the segmentation module were voxelized to downsample to ARA 10 urn resolution while avoiding information loss. Subsequently, the voxelized segmentation maps were warped to ARA 10 urn using a deformation matrix obtained by the registration framework. Middle panels show an example subject overlaid on ARA labels from each group. Lower panels show a 3D rendering of voxelized and warped segmentation maps color-coded based on 6 ARA regions: CTX, Cerebral Cortex; CNU, Cerebral Nuclei; MB, Midbrain; HB, Hindbrain; IB, Interbrain; and CB, Cerebellum. FIG. 6C shows an independent student t-test was applied between c-Fos+ cell density per ARA label to perform a whole-brain ROI-wise statistical test. Upper left panels show the trending ROIs (p < 0.1) while the right panel shows only significant regions (***p < 0.001, **p < 0.01, and *p < 0.05). Lower panels demonstrate the projecting p values on ARA labels. FIG. 6D shows an illustration of the voxelized and warped segmentation maps in each group being averaged and then subtracted toobtain group-wise heatmaps of c-Fos+ neuronal density, identifying neural activity hotspots. FIG. 6E shows the result of the cluster-wise analysis. Upper panel showcases significant clusters inside the midline group of the dorsal thalamus (MTN) including one cluster close to the Xi region between PVH and above the third ventricle (left panel), and multiple clusters in the dorsal subregion of the paraventricular nucleus of the thalamus (PVT; right panel). Middle panels demonstrate the Pearson correlation (only significant correlations are shown, p < 0.01, bottom left triangle) and the Euclidean distance (top right triangle) between the centroid of 22 significant clusters of neuronal ensembles differentially activated after cold stress. The clusters are ranked according to the strength of activations (F- statistic), cluster-wise p-value, and volume. The arrow points to a significant correlation between a cluster mostly located in the nucleus of reunions (RE) and another one in the nucleus accumbens (ACB). Lower panels show 3D views of connectivity / projection maps obtained by anterograde viral vector (AAV) tracing from an experiment (I D= 184158290) in the Allen connectivity atlas, showcasing structural connectivity between RE in MTN and ACB. Left to right panels show MTN and ACB brain regions in 3D, fibers originating from RE and projecting into ACB, maximum intensity projection of signal intensity (highlighting AAV injection site and ACB regions where fibers are mostly traveling into), and a sagittal view of the atlas (cursor is pointing to the ACB region), respectively.

[0078] In another example experiment, identifying locomotive behaviour was examined. Locomotive behaviour is a complex process that involves coordinated neuronal activity in a slew of different brain areas. Identifying activated laminar and layer-specific ensembles of neurons across the entire brain that underlie locomotion represents an ongoing substantial challenge in neuroscience. Using I EG markers (c-Fos), iDISCO+ technique, and CLARITY-optimized light sheet microscope, were used to map the brain regions that contribute to walking in rodents, and spinal cord-projecting neurons underlying the recovery of walking after incomplete spinal cord injury (SCI). Using an ROI-based strategy, elevated levels of c-Fos expression in brain regions were determined that are known to contribute to the production of walking, including the cuneiform nucleus, pedunculopontine nucleus (CUN), primary motor cortex (MOp) and medullary reticular formation. Moreover, using adeno-associated virus tracing and photostimulation, it was demonstrated that a subpopulation of neurons known as the Vglut2 in the lateral hypothalamus (LH; Lateral Hypothalamic Area in ARA) plays a crucial role in walking recovery after chronic lateral hemisection spinal cord injuries (SCI). LH is a specialized brain region primarily recognized for its pivotal role in regulating arousal, feeding, and motivated behaviors, yet it has also been shown to be involved in the orchestration of movement. The ROI approach required performingextensive manual parameter tuning and post-processing using Arivis Vision4D and Fiji plugins for each subject to segment and quantify c-Fos positive cells; using the highest threshold for detection (i.e. the most stringent criteria to minimize false positives), potentially underestimating true cell densities (higher rate of false negatives). Using an ROI-based analysis and a priori hypotheses regarding key areas that are activated in response to walking, it narrows down the c- Fos analysis to known ROIs (contributing to locomotion) and large, homogeneous regional effects.

[0079] In the example experiments investigating locomotion, segmentation maps obtained from the DL models were checked by inspection of voxel-wise density and uncertainty maps, along with the raw data. These segmentation maps were then voxelized and warped (using deformation fields obtained from registration) to the 10pm ARA. By averaging voxelized and warped segmentation maps, the experiments generated a group-wise neuronal density heatmap, demonstrating activity hotspots during walking. The experiments computed c-Fos density per ARA atlas label and performed a whole-brain ROI-based comparison between the two groups. As expected, in the walking group, c-Fos density showed significant increases in the primary motor areas (MOp, and secondary motor area, MOs) in the walking vs homecage group (p < 0.05). Notably, the LH and the midbrain reticular nucleus were among the regions that trended toward significance (p < 0.1), reflecting a localized increase in c-Fos activation that was missed by the ROI-based whole-brain screen, in the uninjured walking vs. homecage (rest) dataset. In addition, ROI-based comparison analysis was repeated at depth 6 instead of the whole brain by merging ARA regions using the hierarchical structure (ontology) of the atlas. Although this analysis identified major areas (at the top of the atlas hierarchy) that contain significantly elevated densities of c-Fos+ cells such as the somatomotor areas (MO, p < 0.05), it missed regions with local and laminar neuronal activations, e.g. the primary somatosensory area (SSp, p > 0.05). To test whether the system 100 is able to detect sub-regional changes, the experiments employed the cluster-wise threshold-free cluster enhancement (TFCE) permutation approach on the voxelized-warped segmentation maps. The outputs of the system 100 highlighted layer-specific areas of c-Fos activation in primary and secondary motor areas and identified several laminar neuronal ensembles differentially activated in walking mice. Interestingly, the example experiments located clusters that are confined to only one layer of the somatomotor areas such as MOp Layer 6a, where thalamocortical projections from the motor thalamus were observed using anterograde viral tracing, as well as clusters that spanned multiple layers within the MOp and Retrosplenial areas. Notably, in contrast to ROI-based analysis, the system 100 detected twodistinct clusters of c-Fos positive neuronal activation in LH, demonstrating their role in movement. To validate the significant cluster, the experiments warped back the registered LH label from the atlas to the native space (using the deformation matrix obtained via registration). Through this analysis, it was determined that higher c-Fos+ cell densities in the walking group vs. homecage, confirming the role of the LH in movement. These findings demonstrate that the system 100 can robustly identify local and laminar brain-wide activity of specific neuronal ensembles in response to experimental manipulation.

[0080] FIGS. 7A to 7F illustrate the locomotion example experiments for brain-wide identification of local neuronal activity changes underlying walking. FIG. 7A shows an overview of experimental design to analyze c-Fos+ cell distribution in whole-brain LSFM data during walking (n = 3 / group). FIG. 7B shows automated segmentation of c-Fos+ cell distribution using the segmentation module. Panels show a 3D Ml P rendering of raw data from the walking group. FIG. 7C shows segmentation maps that were voxelized to the ARA 10 urn resolution while avoiding information loss. Subsequently, the voxelized segmentation maps were warped to ARA 10 urn using a deformation matrix obtained by the registration framework. Left panels show an example downsampled subject overlaid on ARA labels after registration from each group. Right panels show a 3D rendering of voxelized and warped segmentation maps color-coded based on 6 ARA regions: CTX, Cerebral Cortex; CNU, Cerebral Nuclei; MB, Midbrain; HB, Hindbrain; IB, Interbrain; and CB, Cerebellum. FIG. 7D shows an identification of neural activity hotspots by group-wise heatmaps of neuronal density by subtracting the average of the voxelized and warped segmentation maps in each group. Panels show two different coronal views as an example. FIG. 7E shows results of cluster-wise threshold-free cluster enhancement permutation analysis, using a group-wise ANOVA. The panels demonstrate the resulting p-value map representing the clusters showing significant differences between groups and corresponding to the coronal sections in FIG. 7D. Zoomed views show two significant clusters in MOp (left panel) and Retrosplenial area (right panel). FIG. 7F shows that a Lateral Hypothalamic Area (LHA) label was warped back into native space using the deformation matrix obtained by registration. Left and right panels show two example subjects from the walking and homecage (non-walking) groups respectively with a zoomed version of LHA, showcasing higher c-Fos+ activity in the walking vs. homecage condition.

[0081] As illustrated in the example experiments, embodiments of the present disclosure provide an end-to-end pipeline for mapping neuronal ensembles at the sub-regional and laminar levels in tera-voxel scale LSFM data of brains. Embodiments of the present disclosure can leverage ViT-based 3D DL segmentation models, trained on datasets that can be, for example, five times larger than other deep learning approaches for cell segmentation in LSFM. The example experiments validated the segmentation approach using out-of-distribution unseen datasets, showcasing its generalizability. The advanced TFCE permutation-based statistical algorithms for cluster-wise analyses, along with integration with registration approaches, enable automated and sensitive detection of changes in neuronal subpopulations. To demonstrate the versatility and impact, the present embodiments were applied to two distinct neurobiological contexts. Firstly, the system 100 revealed local clusters of neuronal ensembles orchestrating cold-stressed food-seeking behavior. The segmentation and cluster-wise analysis facilitated identifying and validating the role of the Xiphoid region, a small midline thalamic nucleus lacking a predefined atlas boundary in ARA, as a key region differentially activated during food-seeking following prolonged cold stress. Secondly, using the system 100, several sub-regional and laminar neuronal ensembles were discovered to activate during movement, including clusters confined to the secondary motor area, layer 2 / 3, LH regions, or those spanning different layers of MOp. This finding propels the understanding of locomotion as a complex process involving coordinated neuronal activity across various areas.

[0082] In the example experiments, the large training dataset and cutting-edge architectures, in addition to the extensive data augmentation and optimization, enabled the DL models to learn variations in signal-to-noise ratio (SNR) / contrast-to-noise ratio (CNR) levels and intensity profiles of labeled cells, and imaging artifacts across different ROIs and subjects. The system 200 was determined to surpass other image processing-based techniques that require extensive parameter tuning. Unlike other pipelines, the segmentation approach utilizes ensembling via the MC dropout approach to provide a quantitative estimate of model confidence (uncertainty), and increase generalizability; as demonstrated by the system’s 100 ability to segment neurons in out- of-distribution whole-brain datasets and generate robust results across different ROIs. The experiments on the unseen dataset indicated that, on average, the model was up to 50% more precise compared to state-of-the-art segmentation and detection algorithms. Due to the modular nature of the pipeline, the system 100 is readily extendable to different imaging modalities, such as multi-photon microscopy data and 3D histology stacks, via fine-tuning the models using modules that are readily available within the pipeline. Notably, in addition to mapping neuronal activity, the system 100 could be employed to study connectivity via quantifying neuronal soma in upstream regions using retro-grade adeno-associated virus (AAV) viruses.

[0083] The clusters outputted by the system 100 comprehensively characterize significant activity ‘hotspots’, brain-wide, while providing better sensitivity and statistical power compared to those of ROI-based and voxel-wise analyses. The cluster-wise permutation analysis advantageously leverages neighborhood information, mitigating the need for data-specific expertise and minimizing bias in statistical analysis.

[0084] Also advantageously, the system 100 is able to account for covariates, such as sex, age, or weight, or perform mixed effects modeling at the cluster or voxel level. The system 100 can account for covariates by integrating the covariates into the statistical determinations. While the present disclosure generally describes use of ANOVA, it is understood that the system 100 can instead use a linear regression model or use ANCOVA (Analysis of Covariance) to identify statistical differences between the means of two or more independent groups after controlling for one or more explanatory variables (covariates). For examples, covariates can be included in the voxel-wise statistical determinations, and thus, used to refine cluster definitions by considering their effects.

[0085] In the example experiments, neuronal activity associated with cold-induced food-seeking behavior was examined. When studying food-seeking after cold stress, the results aligned with ROI-based semi-automated studies of c-Fos activations. However, the system 100 was able to uncover localized changes in neural activity that likely would have gone unnoticed in ROI-based analyses or required extensive manual inspection to detect. ROI-based analyses may overlook subtle changes as effects, and hence significance, are estimated by averaging over the entire region. In contrast, the system 100 provides a more granular understanding of neural activity patterns, offering enhanced specificity to detect nuanced changes in brain activity over very small areas. The system 100 automatically identified a small neural ensemble associated with food seeking in the Xiphoid nucleus; the activation that was erroneously associated with the entire RE when using ROI-based methods. Employing cluster-wise connectivity analysis, there were significant correlations (putative connectivity) between clusters within the midline group of the dorsal thalamus and the nucleus accumbens, a structurally connected area to Xi, validated by anterograde viral tracing.

[0086] Gaining circuit-level insights based on specific neuronal populations is important for unraveling the organizational principles of the motor system and comprehending its functionality. The example experiments on locomotion-elicited neuronal activity identified two distinct local clusters of activations inside the lateral hypothalamus (LH). The increased activity was restrictedto a sub-region localized within the LH, which likely contributed to the failure of ROI-based analysis to detect LH activation or remained un-noticed in studies utilizing a more limited approach of narrowing down the c-Fos analysis to a few key ROIs associated with locomotion or combining ROIs into larger (and fewer) areas. Furthermore, the system 100 successfully identified several laminar clusters of activations inside the primary motor cortex (MOp), with some clusters being confined to only a single layer, while others spanned multiple layers. This distinctive capability means that the system 100 surpasses the limitations of ROI-based methods and facilitates a comprehensive and unbiased exploration of the intricate organization of the motor system or other brain networks. Thus, the system 100 offers a more nuanced understanding of layer-specific neural dynamics coordinating behavior, especially in regions without existing atlas sub- parcellations. In this way, the system 100 provides a quantitative approach to map locally activated regions in an atlas-agnostic fashion, which is particularly beneficial in experiments exploring brain areas that lack standard or high-resolution digital atlas parcellations. The resulting maps and information about significant clusters can serve as a basis for creating and refining high-fidelity functional or structural atlases. Additionally, leveraging the cluster-wise statistical module and correlation analysis, the system 100 facilitates the extraction of associations and patterns among significant clusters of activations, potentially revealing long-range connectivity or functional coherence.

[0087] In the example experiments, the training data consisted of two LSFM cohorts of TRAP2- Ai9 mice (total N=18): ten animals at the age of three to four months with data acquired from the whole brain; and eight animals at the age of two months with data acquired from the left hemisphere. These animals had an immediate early gene (c-Fos) promoter, contributing to the activity-dependent expression of inducible recombinase Cre-ERT2. The expression of fluorescent protein tdTomato was turned on upon 20mg / kg 4-hydroxytamoxifen (4-TM) injection. The following splits for training and evaluation were used: the training set consisted of ten animals, with three and five animals used for the validation and test sets, respectively.

[0088] In the example experiments, the unseen dataset used whole-brain LSFM data of a model of spinal cord injury (N=3 mice) to validate the segmentation models. Briefly, mice were induced under anesthesia with a mixture of isoflurane and 02. The mice were then placed on a heating pad set to 37°C and maintained on 1-3% isoflurane and 02. A midline skin incision was made and the T10 lamina identified. A T10 laminectomy was performed followed by a left lateral hemisection using micro scissors or a microscalpel as described previously. Muscle closure was performed with 6-0 Vicryl followed by skin closure with 6-0 Ethilon. Post-operatively, mice wereplaced on a heating pad and given subcutaneous fluids as needed. Pain control for 48 hours post- operatively was also provided via subcutaneous daily administration of Rimadyl (5 mg / kg). Bladders were expressed twice daily until spontaneous recovery of bladder function. All hemisection lesions were histology confirmed to be adequate. All three mice were adult female C57BL / 6 mice (>8 weeks at the start of the experiment, 15-30 g body weight).

[0089] In the example experiments, the walking dataset used two groups of adult female C57BL / 6 mice (>8 weeks at the start of the experiment, 15-30 g body weight, n=3 / group). Briefly, mice were trained to run on a treadmill quadrupedally 5 days a week for 2 weeks prior to perfusion. To elicit c-Fos expression, mice ran on the treadmill for 45 minutes at a speed of 9 cm / s. Mice were then perfused 1 hour after (“walking group”). Mice in their home cages were perfused for the “nonwalking group” analysis.

[0090] In the example experiments, fortissue clearing and immunolabeling of the training dataset, two weeks after the 4-TM injections, mice were perfused; their brains were collected and underwent overnight fixation in a 4% paraformaldehyde (PFA) solution. The fixed brain specimens were treated with the SHIELD protocol (Stabilization Under Harsh Conditions via Intramolecular Epoxide Linkages to Prevent Degradation) to maintain the integrity of protein antigenicity. Thereafter, an active clearing procedure was used for tissue CLARITY. The samples were index- matched by adding them to Easylndex medium before imaging them with a light sheet microscope.

[0091] In the example experiments, for tissue clearing and immunolabeling of the unseen and walking datasets, adult mice were anesthetized with intraperitoneal pentobarbital (150 mg / kg) which was followed by intracardiac perfusion of 1 x PBS then 4% PFA in PBS. The brain was dissected, and the sample was post-fixed in 4% PFA overnight at 4°C. Brains then underwent processing with iDISCO+. Briefly, samples underwent methanol pretreatment by dehydrating with a methanol / H2O series, each for 1 hour, as follows: 20%, 40%, 60%, 80%, and 100%. The samples were then washed with 100% methanol for 1 hour then chilled at 4°C followed by overnight incubation in 66% dicholoromethane / 33% methanol at room temperature. This was followed by two washes in 100% methanol at room temperature and then bleached in chilled fresh 5% H2O2 in methanol overnight at 4°C. The samples were rehydrated with a methanol / H2O series as follows: 80%, 60%, 40%, 20%, and PBS each for 1 hour at room temperature. The samples then underwent washing for 1 hour x 2 at room temperature in PTx.2 buffer and then incubated in permeabilization solution for two days at 37°C. Samples were then incubated in blockingsolution (42 mL PTx.2, 3 mL of normal donkey serum, 5 mL of DMSO for a total stock volume of 50 mL) for 2 days at 37°C with shaking. This was followed by incubation in primary antibody solution consisting of PTwH, 5% DMSO, 3% normal donkey serum, and c-Fos (rabbit anti-cFos, 1 :2000, Synaptic Systems, #226 003) for 7 days at 37°C with shaking. After, the samples were washed in PTwH for 24 hours and then incubated in a secondary antibody solution consisting of PTwH, 3% for 7 days at 37°C with shaking. The samples were then washed in PTwH for 24 hours followed by tissue clearing. Final clearing was performed as described. Briefly, samples were dehydrated in a methanol / H2O series as follows: 20%, 40%, 60%, 80%, 100% x 2 each for one hour at room temperature. This was followed by a 3-hour incubation in 66% dichloromethane / 33% methanol at room temperature then incubation in 100% dichloromethane for 15 minutes x 2. The samples were then incubated in dibenzyl ether for at least 24 hours before imaging.

[0092] In the example experiments, light sheet fluorescence microscopy (LSFM) imaging was performed on the training dataset by index-matching the samples with Easylndex medium and then imaged on a light sheet microscope (SmartSPIM) with a 4x objective lens. For eight subjects, data was acquired from a hemisphere while for the rest of the subjects, data were acquired from the whole brain. Nominal lateral spatial resolution was 3.5 pm, and the step size was set to 4 pm in the z-direction with 561 nm excitation and an emission filter of 600 / 52 nm.

[0093] In the example experiments, light sheet fluorescence microscopy (LSFM) imaging was performed on the unseen and walking datasets by imaging of whole brains using a COLM (CLARITY-optimized light sheet microscope) with a pixel resolution of 1.4x1.4 pm in the x, y dimensions and a z-step of 5 pm, respectively, using the 4x / 0.28 (d) objective which is suitable for resolving cell nuclei labeled by c-Fos. Using a custom-made quartz cuvette filled with dibenzyl ether, whole brains were imaged. One channel with autofluorescence (“auto” channel; 488 nm) to demonstrate anatomy was imaged and the other channel demonstrating cFos labelling (“cell” channel; 647 nm) was imaged. All raw images were acquired as 16-bit TIFF files and were stitched together.

[0094] In the example experiments, ground truth (GT) labels of neuronal somas comprised three stages: segmentation, llastik pixel classification, and post-processing. To create GT labels, the segmentation workflow was performed that incorporated image processing tools implemented as FIJI / lmageJ macros. The workflow includes a 3D Watershed marker-controlled algorithm and a post-processing 3D shape filter to omit false positives (FP). This resulted in binary GT labels across the entire dataset with a low FP rate and a relatively higher false negative (FN) rate.

[0095] To improve the GT labels, the example experiments used llastik, which performs pixel classification using image filters (as input features) and a random forest (RF) algorithm (as a classifier). Filters include pixel color and intensity descriptors, edginess, and texture in 3D and at different scales. The RF combines hundreds of decision trees and trains each one on a slightly different set of features. The final predictions of the RF are made by averaging the predictions of each tree. The example experiments used 37 3D filters and an RF classifier with 100 trees. To train the RF classifier, the segmentation outputs were imported as “silver” input annotations (initialization) to llastik (i.e. , in lieu of manual annotations). To address the prohibitive speed and memory requirements of the RF algorithm, it was trained on image patches (5123voxels). For each brain, three 5123patches from different depths were randomly selected. The RF model was trained several times by providing feedback (i.e., by correcting the results with expert annotation and modifying the labels to reach the best results). The llaskik output did not correctly detect the boundaries of neuronal soma and frequently overestimated their spatial extent.

[0096] In the example experiments, to further reduce the number of FP voxels in the GT labels, a 3D shape filter was applied using the Imaged shape filter plugin on the llastik-generated labels. To solve the volume overestimation problem, a 3D erosion filter was applied with a sphere-like kernel (radius of 1 voxel) on the GT labels. The final whole-brain GT labels were visually quality- checked in three 5123patches per brain by two raters. Specifically, the raters were asked to randomly select and validate one 5123patch in the cerebrum, brain stem, and cerebellum.

[0097] In the example experiments, the DL models implemented by the segmentation module 104 included 3D ViT-based (LINETR) and CNN-based (ll-Net) architectures.

[0098] LINETR has powerful representation learning capabilities and generates more accurate dense segmentation masks than do other CNN architectures (such as Mask R-CNN), thanks to the preservation of the spatial information. However, fully convolutional models are limited in their ability to learn long-range dependencies; resulting in potentially sub-optimal segmentation of objects in large volumes, including neuronal cell bodies with different shapes and size (e.g., in different regions of the brain). To address this issue, a ViT architecture was used in which the encoder path of a ll-Net is replaced by a transformer to learn contextual information from the embedded input patches. In this way, the UNETR-based segmentation models had residual blocks and dropout layers to improve the robustness and generalizability of previous pipelines. Specifically, in the LINETR, the encoder was replaced with a stack of 12 transformer blocks, operating on a 1 D sequence (163= 4096) embedding of the input (one channel 3D image patch,x e / ?96x96x96xi Subsequently, a linear layer projected the vectors into a higher dimensional embedding space (K = 768), which remained constant throughout the transformer layers. A 1 D learnable positional embedding was added to the projected patch sequence to preserve the spatial information. This embedded patch sequence with the dimension of N x K (N = ^ x ^ x ^ = 216 is the sequence length) was passed as the input to the stack of transformer blocks. Each transformer block consisted of multi-head attention and multi-layer perceptron (MLP) layers. Each multi-head attention module consisted of multiple self-attention heads, where each self-attention block learned the mapping in the patch sequence in parallel. Learned sequence representations were learned at 4 different depths of the transformer stack and reshaped eachback to — x — x — x K. These vectors then went through five different encoders with consecutive16 16 16aconvolutional layers to achieve supervision at different depths (voxel sizes).

[0099] In the UNETR of the example experiments, the encoder was connected to a decoder via skip connections at multiple resolutions to predict the segmentation outputs. Dropout layers were included in all blocks and residual units were included in convolution blocks in the encoder as a regularization strategy to avoid overfitting and vanishing gradient problems. For a given 3D input cube, the segmentation models generated a 3D volume containing voxel-wise probabilities of neuronal cell bodies (0 < P < 1). Prediction maps were binarized (using a threshold of 0.5) to generate a map representing whether a voxel belongs to a neuron or not.

[0100] With respect to the U-Net of the example embodiments, it consists of a contracting (encoder) and expanding (decoder) paths. The encoder is based on 3D convolution and pooling operators. It takes an image patch as input and generates feature maps at various scales, creating a multi-level, multi-resolution feature representation. Meanwhile, the decoder with up-convolution operators leverages these feature representations to classify all pixels at the original image resolution. The decoder assembles the segmentation, starting with low-resolution feature maps that capture large-scale structures and gradually refining the output to include fine-scale details. In the U-Net architecture, the standard building blocks have been replaced with a residual block in the example experiments. In addition, parametric rectifying linear units have been deployed to provide different parameterized nonlinearity for different layers. The 2D operations have been included that include convolution and max-pooling layers into 3D and use batch normalization instead of instance normalization to achieve a stable distribution of activation values throughout the training and to accelerate the training. A dropout layer has been added between all convolution blocks in the architecture as a regularizer to avoid overfitting. Except for the first residual units, allthe convolutions, and transpose convolutions had a stride of 2 for downsampling, and upsampling, respectively. The first residual unit used a stride of 1 which has been shown to increase the performance by not immediately downsampling the input image patch.

[0101] With respect to the loss function used in the example embodiments, Dice coefficient is a metric that measures the similarity between two labels. The class average Dice can be computed as:Where N is the number of classes, K is the number of voxels, and GTk nand Pk ndenote the one-hot representation of ground truth and the output prediction maps for class n at voxel k, respectively.

[0102] Cross entropy measures the difference between two probability distributions over the same sets of underlying events and was computed as:

[0103] The example experiments used an equally weighted Dice-Cross Entropy loss, which is a combination of Dice loss and Cross Entropy Loss functions. It was computed in a voxel-wise manner as:L(GT, P') = 1 - (D(GT, P) + CE(GT, P))

[0104] With respect to model training used in the example embodiments, a total of 36,480 (963voxels) unique input patches were used, not accounting for data augmentation. To address the class imbalance issue in the dataset (majority of voxels representing background), patches containing only background, or a small number of foreground voxels (<100,000 foreground voxels) were filtered out from the 5123image patches generated by the annotation strategy. Hence, for the LINETR model with an input size of 963, 15,200, 9120, and 12,160 unique input patches (not accounting for data augmentation) were used for training, validation, and test (N=30,400 patches with augmentation for training alone). While for the ll-Net model with an input size of 1283, 7600, 3840, and 5120 patches were used respectively.

[0105] For hyperparameter tuning in the example embodiments, a Bayesian optimization approach was used. The following hyperparameters were optimized during training: input image size, encoder and decoder’s depth, kernel size, learning rate, batch size, size of kernels, and loss function. For the LINETR architecture, which in total had 92.8M parameters, the best model based on Dice score performance on the validation set had the following parameters: 12 attention heads, a feature size of 16, an input patch size of 963voxels, and a batch size of 24. The ll-Net architecture, which in total had 4.8M parameters, had the following parameters: a 5-layer encoder with an initial channel size of 16 and a kernel size of 3x3x3, input patch size of 1283voxels, and batch size of 27. The LINETR model was trained for 700 epochs, while the ll-Net model was trained for 580 epochs. Early stopping was set to 50 epochs where performance (Dice score) on the validation dataset did not improve to avoid overfitting. The Adam optimizer was used with an initial learning rate of 0.0001.

[0106] To increase the generalizability of the DL model and model distribution shifts in LSFM data, in the example experiments, the training set images were randomly augmented in real time at every epoch;, namely: affine transformations, contrast adjustment, histogram shift, random axis flipping, and different noise distributions such as salt and pepper, and Gaussian. These transforms were selected as they are representative of distortions that occur during LSFM imaging. Briefly, a range (0, 1) scaling was applied to each 5123image patch based on the intensity distribution of the patch. The intensity of each 5123image patch was scaled from [0.05, 99.5] percentile to [0, 1], where 0.05 and 99.95 are the intensity values at the corresponding percentiles of the image patch. Subsequently, each data augmentation transform was randomly applied (with a probability of p = 0.5) to the 5123image patch at each epoch. The parameters of each data augmentation were also randomly selected at each epoch from the predetermined range of values.

[0107] In the example experiments, a voxel-wise uncertainty map was generated. Epistemic uncertainty is rooted in a lack of knowledge about the model parameters and structure, rather than stemming from the inherent variability in the observed data. This type of uncertainty is often referred to as model uncertainty. To estimate the models’ uncertainty and confidence in predictions, a Monte Carlo (MC) dropout approach was used. During training, dropout layers (with a probability of p = 0.2) were utilized as a regularization technique. It has been shown that turning on dropout layers (randomly switching neurons off) in inference mode can be interpreted as a Bayesian approximation of the Gaussian process. At test time, when the dropout layers are turned on, each forward pass yields a stochastically different prediction as a sample from theapproximate parametric posterior distribution (P(y |X)). This technique has been shown to provide useful insights into the model's uncertainty by computing the variance of numerous predictions (T = {yi,y2, ...,yjv})- Voxel-wise uncertainty (variance) can thus be defined as:Voxel-wise uncertaintyThe resulting voxel-wise uncertainty map provides a measure of how much predictions of slightly different models vary on the same data, which can be particularly useful in identifying regions of high uncertainty.

[0108] For the ensemble of ensembles approach used in the example experiments, each model’s output was obtained by averaging 50 models’ probability maps (T = {y1,y2, . .. ,y50}) using the MC dropout technique:Subsequently, to combine both models’ prediction maps, a final mapping of neurons was generated using an ensemble of both models (ensemble of ensembles):

[0109] In the example experiments, the models were evaluated using evaluation metrics. To evaluate the performance of the segmentation models, several volume and shape-based metrics were used including the Dice similarity coefficient (DSC), recall, precision, F1 -score, and %95 Hausdorff distance. The metrics were derived from the resulting confusion matrix and the associated true positive (TP), false positive (FP), and false-negative (FN) values. Each neuronal soma was defined as a 3D connected component. Given this definition, TP was defined as the number of correctly detected neurons after comparing the ground truth GT to the prediction P. Sensitivity or Recall measures the proportion of TP to the number of individual neurons delineated in GT and was defined as:TPRecall =TP+FNPrecision measures the proportion of TP against all positive predictions and was defined as:TPPrecision =TP+FPWhile Recall is useful to gauge the number of FN pixels in an image, Precision is useful to evaluate the number of FP pixels in a prediction. The F1-score combines Precision and Recall and is often used to measure the overall performance of a model. The F1 -score measures the number of wrongly detected neurons in P:2xPrecisionxRecall F1 -score = -Precision + Recall

[0110] The DSC metric measures the number of elements common to GT and P datasets divided by the sum of the elements in each dataset. The DSC was defined as:

[0111] Hausdorff distance metric is a mathematical measurement on the “closeness” of two sets of points that are subsets of a metric space. It is the greatest of all the distances from a point in one set to the closest point in the other set. The example experiments used the 95 percentile of the Hausdorff Distance rather than the maximum results to provide a more robust and representative measure of the segmentation's performance in image analysis tasks. Given two sets of points X and Y, the Hausdorff distance between these two sets was defined as:HD(X, Y) = max {supxEXinfyEY^x> y)> supyEYinfxEX^x> y)}

[0112] For the comparison against llastik in the example experiments, the pixel classification module in llastik was used and trained a random forest (RF) classifier using all training subjects (18 whole-brain LSFM images). 37 3D filters and an RF classifier with 100 trees was used. To train the RF model, the experiments randomly selected three 5123voxels image patches from each subject and applied the same scale-intensity transform as used in the DL training approach to provide a fair comparison. The trained RF classifier was used to generate segmentation maps for all test and unseen datasets. To quantitively evaluate the system’s 100 robustness across different regions of the brain, the test set was registered to ARA 10 pm, warped ARA labels back to their native space, and integrated warped labels with segmentation maps.

[0113] For the comparison against Cellfinder in the example experiments, Cellfinder was used to generate whole-brain detection maps of neuronal cell bodies. Their pre-trained Resnet model was used to generate the detection maps on both test and unseen datasets. The cellfinder commandwas deployed with -no-registration flag to detect and classify cells to background or neuron. For retraining the Cellfinder, there was random selection of two subjects from the test set and generated around -6100 (-3200 cell and -2900 non-cell) annotated cell candidates for the first subject and -9000 (-4300 cell and -4700 non-cell) for the second subject using the Napari- cellfinder plugin. The training data was then used to retrain the Resnet model by incorporating cellfinder_train function and -continue-training flag, keeping other options as default. Lastly, the best retrained model based on validation error was used to generate whole-brain detection maps. To quantitatively compare against Cellfinder, the segmentation maps generated by the system 100 were transformed into detection maps by finding the centre of mass of each neuron in 3D.

[0114] To synthesize and correlate segmentation results within the ARA space, the example experiments employed a voxelization approach. Voxelization entailed the transformation of high- resolution segmentation outcomes into a 10-pm resolution space while avoiding loss of information. Segmentation volumes underwent convolution with a spherical kernel featuring a -5 pm radius; a dimension that aligns with the downsampling factor and ARA space. Subsequently, within each convolved sphere, the average count of labels (cells or nuclei) was computed, resulting in a voxelized map. This 3D voxelized representation, which was generated using Python's skimage library with parallel computation facilitated by the joblib and multiprocessing libraries, allowed for efficient feature extraction summarized by ARA regions / labels.

[0115] To bring whole-brain tissue-cleared microscopy images, segmentation maps, and a reference atlas (Allen Mouse Brain Atlas, ARA) into spatial correspondence, the example experiments used specialized workflows that were optimized for multi-modal registration of cleared data. The registration workflows include a cascade of image preprocessing techniques such as denoising and intensity correction, as well as an intensity-based alignment. The alignment process consisted of two main steps. In the first step, an initial alignment is carried out using the antsAffinelnitializer tool from Advanced Normalization Tools (ANTs). The second step consisted of an intensity-based multistage b-spline registration algorithm, encompassing a rigid 6 degrees of freedom (DOF), an affine (12 DOF), and a non-rigid (deformable) b-spline symmetric normalization stage. The rigid and affine stages are based on maximizing the mutual information similarity metric between ARA and microscopy data, while the deformable stage used crosscorrelation as the similarity metric. The resulting transformations perform bidirectional warping of images to and from tissue-cleared microscopy native space and ARA space.

[0116] In the example experiments, voxelized segmentation maps were warped to the ARA space (10 pm resolution) using a deformation field obtained through registration. Subsequently, a smoothing Gaussian filter was applied with a sigma of four pixels. Difference heatmaps were then computed by subtracting the average of voxelized and warped segmentation maps in each group.

[0117] In the example experiments, the voxelized and warped segmentation map, at a resolution of 10pm, in addition to the registered labels for each subject were used to extract the density of cells per ARA region, for whole brain as well as labels grouped to a maximum atlas ontology depth of 6. The ARA labels are structured into descending depths going from coarse to finegrained groupings of brain regions with a function that combines labels at a higher depth by their parent labels. The resulting density results were then passed to a function which applied Student's paired t-test (two-sided) per label with alpha value of 0.05 for both whole brain regions and depth 6 regions and both applications. The function creates bar plots that compare the density of cells per ARA label including both significant ROIs and trending regions (p-value < 0.1). Lastly, a function was used to project the resulting p-values on the atlas regions.

[0118] A cluster-wise permutation-based statistical analysis with threshold-free cluster enhancement (TFCE) was used by the system 100 in the example experiments. The statistical pipeline consisted of three main steps. In the first step, a voxel-wise statistical test using a two- way ANOVA was performed between two groups of study. To incorporate the correlation structure of the data and correct for multiple comparisons, a null hypothesis was considered regarding the sizes of clusters in the data instead of focusing on individual voxels. Thus, in the second step, clusters were defined using adjacency structure in the data (connecting each voxel to its neighbors in 3D) and the TFCE technique was applied, which addresses the challenge of selecting a threshold for cluster forming as well as smoothing problems. The example experiments optimized the adjacency structure of the data using a priori knowledge from the group-wise heatmaps to boost the sensitivity of the statistical results in tera-voxel LSFM data, due to the exceedingly large number of voxels (and hence the number of statistical tests). Specifically, the adjacency matrix was masked using a thresholded (90% percentile) then dilated version of the group-difference heatmap, allowing the algorithm to focus only on putative clusters. TFCE transformed a raw statistical image into a map that reflects the strength of local spatial clustering. The new value assigned to each voxel is determined by aggregating the scores of all supporting sections below it. The score of each section is computed by taking the height of the section (raised to a power H) and multiplying it by the extent of the section (raised to a power E):Where h0is typically around zero and / iris statistic value corresponding to voxel v. In practice, this integral is estimated as a sum, using finite step sizes (dh). The exponents of the powers (E and / 7) are free parameters, but fixing these values yields robust outcomes. Increasing H gives more weight to clusters with higher effect sizes while increasing E gives more weight to larger clusters. In the example experiments, E = 2, H = 0.5, h0= 0.0, and a step-size of 5 was selected to mitigate the potential impact of false positives in LSFM data and fixed the values for both applications, walking and food-seeking datasets. Applying TFCE transformed results in a weighted sum of local cluster-like signals, eliminating the need for a fixed threshold to define clusters while keeping the local minima and maxima at the same spot. The size of each cluster was measured by the sum of voxels’ F-values within the cluster. In the last step, a non-parametric permutation test (N = 1000) was applied and new cluster-wise F-statistics were obtained. In each permutation, the same approach was applied to define clusters and compute their statistics, and the largest cluster size was retained. The step-down p value procedure was also used to boost sensitivity while controlling for Family-wise error rate. To test the significance of clusters in our actual data, a null distribution was obtained via permutations. Cluster sizes observed in the actual data were compared with those in the null distribution to calculate p-values. This p-value can be compared to a predetermined threshold (such as an alpha < 0.05 or FDR-corrected p-values) to test for significance. Next, the example experiments applied connected component analysis on the resulting p-value image to summarize significant clusters. Finally, by integrating connected component analysis results, ARA labels, and density heatmaps, the example experiments extracted the center, mean effect size, and volume of each cluster, along with the percent volume within each brain region spanned.

[0119] In the example experiments, the significant clusters identified with the cluster-wise permutation-based statistical analysis were systematically ranked based on their statistical effect size and volume. To streamline the subsequent connectivity analysis, the top 20 clusters from this ordered list was selected for further analysis. Integrating each cluster's location with voxelized and warped segmentation maps, the mean intensity of each cluster per subject was determined in both the treated and control groups. Subsequently, a Pearson correlation test, utilizing the pearsonr function from scipy library, was used between the mean intensity values of each pair of clusters. For each correlation coefficient, an 80% confidence interval was established using the bias-corrected and accelerated bootstrap method with 10,000 iterations, and the lower bound wasselected. Finally, the p-value associated with each correlation test was determined through permutation testing, involving 10,000 permutations. For visualization, significant correlation coefficients (p < 0.01) were plotted alongside a heatmap illustrating the Euclidean distance between each pair of cluster centroids.

[0120] While the example experiments primarily focused on mapping neuronal activity in LSFM data using immediate early genes (lEGs) and c-Fos activity, the system 100 can be used with other fluorescent reporters. Advantageously, the ensemble of ensembles model, as described herein, leveraging CNN-based and ViT-based architectures, can be used to deal with potential false positives on unseen datasets stemming from artifacts or autofluorescence mimics from, for example, lipofusion; increasing segmentation performance and boosting statistical specificity. Advantageously, the system 100 does not require pre-defined cluster-forming thresholds that are required in simple cluster-wise inference approaches because it utilizes free parameters to produce results in which the voxel-wise values represent the amount of cluster-like local spatial support. The example experiments show that this produces robust outcomes for 3D LSFM data analyses, as the same parameters were applicable to both food-seeking and walking experiments.

[0121] Accordingly, the embodiments of the present disclosure provide a streamlined approach for achieving high-fidelity unbiased 3D mapping of local and laminar neuronal ensembles across the entire brain, independent of pre-defined atlas regions. The DL models are able to provide robust and generalizable segmentation of neuronal somas in whole-brain tera-voxel LSFM data. The system 100 implements a robust framework for quantitative identification and discovery of differentially activated localized clusters of ensembles in response to experimental manipulations. The outputs of the system 100 permit examination of neural activity patterns at an unprecedented level of spatial detail, providing valuable insights that may be obscured with traditional ROI-based or voxel-wise analysis. Embodiments of the present disclosure are applicable across a wide range of LSFM datasets and neuroscience paradigms.

[0122] In some cases of the present embodiments, the generated uncertainty maps may be leveraged in a feedback loop to further improve segmentation performance. For example, regions with high uncertainty can be excluded during inference (e.g., by zeroing out high-uncertainty areas) to reduce false positives. Additionally, the uncertainty maps may be incorporated into the training / fine-tuning process by integrating them into the loss function, thereby enabling adapting of the models that emphasizes confident regions while penalizing uncertain predictions. Thisstrategy allows the system to adapt more effectively to variations in input data and enhance overall segmentation accuracy.

[0123] Although the foregoing has been described with reference to certain specific embodiments, various modifications thereto will be apparent to those skilled in the art without departing from the spirit and scope of the invention as outlined in the appended claims.

Claims

Claims1. A method for generating a mapping of local brain cell activity and connectivity, the method comprising: receiving brain imaging data from image captures of a brain; generating a binary segmentation map using a trained deep learning model, the deep learning model taking the brain imaging data as input; performing registration of the brain imaging data to a brain template; voxelizing the segmentation map and warping the voxelized segmentation map using deformations obtained during registration; determining group-wise heatmaps of neuronal density or cellular density from the voxelized and warped segmentation map; generating a map of clusters of between-group differences in neuronal density or cellular density from the group-wise heatmaps of the voxelized and warped segmentation map; and outputting the map of clusters as a representation of local brain cell activity and connectivity.

2. The method of claim 1, wherein the brain imaging data comprises whole-brain imaging data.

3. The method of claim 1, wherein the registration comprises performing one or more of: an initial alignment using an initial affine transformation; a rigid registration that allows only translations and rotations to refine an alignment; an affine registration that incorporates one or more of translations, rotations, scaling, and shearing; and a non-rigid, deformable registration using a b-spline symmetric normalization.

4. The method of claim 1, wherein the map of clusters comprises a p-value map.

5. The method of claim 1 , wherein the map of clusters is generated using at least one of voxel-wise comparisons, cluster definition and enhancement, and permutation testing.

6. The method of claim 1, wherein voxelizing the segmentation map comprises using a convolution filter and warping the segmentation map comprises using a brain reference atlas using the deformations obtained during registration.

7. The method of claim 1, wherein the heatmaps are determined by subtracting an average of the voxelized and warped segmentation map in each group of the group-wise heatmaps.

8. The method of claim 1, wherein the between-group differences are determined by performing cluster-wise threshold-free cluster enhancement permutation analysis between the groups on the group-wise heatmaps.

9. The method of claim 8, wherein the threshold-free cluster enhancement permutation analysis comprises a group-wise two-way analysis of variance (ANOVA) test.

10. The method of claim 1, further comprising outputting auxiliary information about the clusters, the auxiliary information comprising at least one of volume, morphological features, location, and portion of each brain region that is included in each cluster.

11. A system for generating a mapping of local brain cell activity and connectivity, the system comprising one or more processors in communication with a non-transitory storage medium, the non-transitory storage medium comprising instructions for the one or more processors to execute: an input module to receive brain imaging data from image captures of a brain; a segmentation module to generate a binary segmentation map using a trained deep learning model, the deep learning model taking the brain imaging data as input; a registration module to perform registration of the brain imaging data to a brain template, and to voxelize the segmentation map and warp the voxelized segmentation map using deformations obtained during registration; a statistics module to determine group-wise heatmaps of neuronal density or cellular density from the voxelized and warped segmentation map, and to generate a map of clusters of between-group differences in neuronal density or cellular density from the group-wise heatmaps of the voxelized and warped segmentation map; andan output module to output the map of clusters as a representation of local brain cell activity and connectivity.

12. The system of claim 11, wherein the brain imaging data comprises whole-brain imaging data.

13. The system of claim 11, wherein the registration module performs the registration by performing one or more of: an initial alignment using an initial affine transformation; a rigid registration that allows only translations and rotations to refine an alignment; an affine registration that incorporates one or more of translations, rotations, scaling, and shearing; and a non-rigid, deformable registration using a b-spline symmetric normalization.

14. The system of claim 11, wherein the map of clusters comprises a p-value map.

15. The system of claim 11, wherein the map of clusters is generated using at least one of voxel-wise comparisons, cluster definition and enhancement, and permutation testing.

16. The system of claim 11, wherein the registration module voxelizes the segmentation map using a convolution filter and warps the segmentation map comprises using a brain reference atlas using the deformations obtained during registration.

17. The system of claim 11, wherein the statistics module determines the heatmaps by subtracting an average of the voxelized and warped segmentation map in each group of the group-wise heatmaps.

18. The system of claim 11, wherein the statistics module determines the between-group differences by performing cluster-wise threshold-free cluster enhancement permutation analysis between the groups on the group-wise heatmaps.

19. The system of claim 18, wherein the threshold-free cluster enhancement permutation analysis comprises a group-wise two-way analysis of variance (ANOVA) test.

20. The system of claim 11, wherein the output module further outputs auxiliary information about the clusters, the auxiliary information comprising at least one of volume, morphological features, location, and portion of each brain region that is included in each cluster.

Images