Data-driven segmentation of cortical calcium dynamics

Abstract

Demixing signals in transcranial videos of neuronal calcium flux across the cerebral hemispheres is a key step before mapping features of cortical organization. Here we demonstrate that independent component analysis can optimally recover neural signal content in widefield recordings of neuronal cortical calcium dynamics captured at a minimum sampling rate of 1.5×106 pixels per one-hundred millisecond frame for seventeen minutes with a magnification ratio of 1:1. We show that a set of spatial and temporal metrics obtained from the components can be used to build a random forest classifier, which separates neural activity and artifact components automatically at human performance. Using this data, we establish functional segmentation of the mouse cortex to provide a map of ~115 domains per hemisphere, in which extracted time courses maximally represent the underlying signal in each recording. Domain maps revealed substantial regional motifs, with higher order cortical regions presenting large, eccentric domains compared with smaller, more circular ones in primary sensory areas. This workflow of data-driven video decomposition and machine classification of signal sources can greatly enhance high quality mapping of complex cerebral dynamics.

Fig 1. ICA separates calcium data into its underlying signal components.

(A) Recording schematic and fluorescence image of transcranial calcium imaging preparation. Full image with masked regions of interest (ROI, dashed lines) shows regions used in ICA decomposition. (B) Sample video montage of raw video frames after ΔF/F calculation. (C) ICA video decomposition workflow. The ΔF/F movie without the mean is decomposed into a series of statistically independent components that are either neural, artifact, or noise associated (not displayed). Each component has an associated time course from the ICA mixing matrix. Neural components can be rebuilt into a filtered movie (rICA). Alternatively, artifact components can be rebuilt into an artifact movie. Circular panels show higher resolution spatial structure in example in the rightmost components.

Fig 2. ICA decomposition quality is sensitive to recording spatial and temporal resolution.

(A) Distributions for lag-1 autocorrelation (black) and log temporal variance (purple) are displayed for components 1–1200. A dotted line represents the cutoff determined from the distribution in the right panel. In the right panel, a horizontal histogram on the lag-1 autocorrelation with a two-peaked kernel density estimator (KDE) fit reveals a two peaked-histogram, summarized by a barbell line. Group data for each peak, as well as the central cutoff value is summarized by the boxplots on the right (n = 16 videos; from 8 different animals). Examples of neuronal (n; blue), artifact (vascular—v; red/other -o; orange), and noise (gray) ICs are indicated in the variance plot are shown below. Comparative ICs across spatial down sampling are shown in S2 Fig and temporal down sampling show in S3 Fig. (B) 2-peaked KDE fits of horizontal histogram distributions under various spatial down sampling conditions, with barbell summary lines on the right. After spatial resolution decreases beyond 41 μm pixel width (px), this two peak structure collapses, and an x denotes the primary histogram peak. (C) 2-peaked KDE fits of horizontal histogram distributions under four temporal down sampling conditions, with barbell summary lines on the right. sIC global max and min of each spatially (D) and temporally (F) down sampled experiment and their temporal standard deviation of the first 1000 data points (E,G).

Fig 3. Increased duration increases specificity of each IC.

(A) Component stabilization for different length video subsets of six 20- minute video samples. (n = 6 videos from 3 different animals) Individual thin lines show polynomial fit to neural or artifact components under each time condition. Thick lines denote the curve fit of the mean number of components in each category across these six experiments. The group distribution of components at 20 minutes is summarized by the boxplot on the right (n = 16 videos; from 8 different animals). Dashed vertical lines are durations in which examples are shown in D. (B) Temporal standard deviation of the first 1000 data points. (C) sIC global max and min of each varying duration included in each decomposition. (D) Examples of similar ICs between decompositions of varying lengths for both neural components (left sIC with central blue time courses) and artifacts (right sIC with central red time courses).

Fig 4. Class identity cannot be established by any individual extracted feature.

(A) Examples of independent components of neural (n) signal, vascular (v) artifacts, and other (o) artifacts. Components are defined by both the sIC and its temporal fluctuations. Circular windows magnify key portions of the sIC. sIC values represented by colormap from blue to red. Temporal representation is in relative intensity (black time course under the sIC), only 1 minute of the full 20 minutes are shown. (B) A comparison of the number of neural signal (GCaMP: dark blue; controls: light blue) and the artifact components (vascular: red; other: orange) with each animal shown (GCaMP components: N = 12 animals, n = 3851; mGFP components: N = 3, n = 484; aGFP components: N = 3, n = 442; WT components: N = 3, n = 229). (C) Examples of binarization of the sIC. Histogram shows the full distribution of sIC values. The dynamic threshold method to generate binarized masks was used to identify the high sIC signal pixels (yellow) against the gaussian background (blue). Windowed spatial representation shows binarization on the key portions of the sIC. (D) Examples of neural and artifact wavelet analysis shown in the power signal-to-noise ratio (PNR) plots. 95% red-noise cutoff was used to create signal to noise ratio (black dashed lines). (E) Histograms of example spatial metrics derived from GCaMP sIC values, (F) morphometrics from the shape of the binarized primary region, (G) temporal metrics derived from relative temporal intensities, (H) frequency metrics derived from the PNR.

Fig 5. Spatial thresholding and frequency data reliably produce neural metrics.

(A) Individual experiment preparation with corresponding spatial footprints by class of component: GCaMP neural (dark blue), control neural (light blue), vascular (red), other (orange). (B) All model experiments (N = 12) with corresponding centroid location of each of the class of metrics. Histograms show the resulting average distribution of spatial location across the field of view (error bars are standard deviation between experiments). (C) Individual experiment (same as A), where components are sorted by temporal variance. PNR mapped to each component and organized between the classification of components. (D) Main frequencies seen in each component class between each experimental condition. Dotted lines represent the mean dominant frequency within each animal, where the gray around the mean corresponds to the standard deviation of that animal. The color line corresponds to the grand mean between all experiments. (E) Relative position-based variance of the types of components between experiments and transgenic model, shown as the average and standard deviation between experiments. (F) The percent of components that had footprints and frequency data that was above the noise cut-off, separated by component type and experimental condition.

Fig 6. Spatial and morphological metrics are most important to classify components.

(A) Correlation and t-statistic between artifact and neural components for each feature (N = 7, n = 2190). Spatial (circles), morphological (triangle), temporal (diamond), and frequency (square) metrics plotted. Cut off values that helped in the selection process are dotted lines, rejected values in gray. Closed points are components that meet requirements. Relative importance metric from the Random Forest classifier plotted against each metric by their respective classes. Selected metrics shown in the list within each type of feature, sorted by greatest t-statistics magnitude. (B) The dataset was parsed into ML modeling dataset (N = 7, n = 2190) that was used to establish the machine learning pipeline and a novel dataset (N = 5, n = 1661) of full experiments that will not influence the classifier. Modeling data was stratified 70/30 split based on classification. 1000 iterations of training the machine learning classifier on selected metrics and validating the machine classification with human classifications. (C) Performance of the ML training, using subsets of the ML modeling dataset. 1000 iterations resulted in accuracy, precision and recall boxplots. (D) 1000 iterations of training on the full ML building dataset was performed and the novel dataset was assessed on its performance. (F) SVD projection of metric data with human classification mapping (top) and the confidence of the ML classifier (bottom). (E) Performance of the classifier on each of the novel datasets, animals plotted separately showing distribution of the 1000 different trained classifiers. (G) Approximate location of false negatives and positives from novel datasets.

Fig 7. Domain maps represent the features of the ICs.

(A) Schematic of domain map creation. A maximum projection is taken through each blurred neural component to form a domain map. Mean time courses are extracted from rebuilt filtered movies within a defined domain. (B) Left scatter plot compares the number of neuronal ICs and the resulting number of domains and the dashed line is the identity line. The circled dot indicates the examples used in this figure. Right bar plot shows the fraction of ICs that have major contributions to the number of domains. (C) Examples of ICs that contributed two domains, one domain or had no contribution. The location of the maxima (black circle) of each IC was found and a point correlation map of the rebuilt full resolution filtered movie was produced (green and pink maps). Pixel-wise scatter plot (right) shows the relationship between the spatial IC value and the correlation map value. Dotted line is the threshold value described in Fig 4. Red line indicates the median correlation value of the correlation map that resides in the thresholded IC. (D) IC-Domain correlation with domain map created with filtered mean plotted with respect to the size of each domain. Dark blue corresponds to domains where each IC made more than one domain. Bright blue corresponds to those that made only one. (E) Median correlation value based on point-correlation analysis within each domain (top). Difference in correlation between the center domain with all its adjacent neighbors (gray scale value corresponds to the mean of the median difference value found in each immediately adjacent domain). (F) Example time courses of domain neighborhood. Left domain map identifies the location of each neighborhood with each corresponding time series. Example ICs and point correlation data shown in S10 Fig.

Fig 8. Time series extracted from domain maps outperform time series generated from other down sampling methods.

(A) Example of a mosaic movie frame rebuilt with respect to each down sampling technique. The non-down sampled filtered movie is represented on the left with subsequent down sampling based on domain, grid or Voronoi maps. (B) Percent total signal of the filtered video represented by extracted time courses. Percent of overall video signal captured in domain maps was calculated for each animal (green circle; N = 8), and compared to signal content from a domain map generated from a separate video from the same animal (green triangle). Percent total signal represented by time courses extracted from grid (blue square) or randomly generated (blue diamond) maps were compared as controls. In the right panel, the percent signal relative to the domain map percent signal was summarized in a box plot. (C) Variation between time courses extracted with each map method was then quantified as a sum signal variation for each experiment. In the right panel, the sum signal variation for each comparison map relative to the optimized domain map sum signal variation was summarized in a box plot.

Fig 9. Domain maps are created from ICA components and are unique to each recording, but highly similar among individual animals.

(A) Domain size (circle size) and minor axis (green-purple colorbar) plotted at each position across all hemispheres independently (left). Highest eccentric domains (top quartile) plotted showing the location and direction of its major axis (right; n hemispheres = 14, 7 mice) (B) The Allen Brain atlas map [36] is additionally used for anatomical reference. (C) The final manually assigned region, with associated labels. (D) Domain minor axis and eccentricity by region (grey horizontal line denotes top quantile used in A). 7 mice). (E) Example overlay of one domain map on another from the same animal. Individual domain or region overlap is calculated using the Jaccard index (intersect / union). Population analysis of the Jaccard index for domain (F) and region (G) overlap comparisons. Maps are generated from a different recording on the same animal, a littermate, a non-littermate, or a randomly generated voronoi map. Significance is calculated using a two-way ANOVA, followed by post-hoc t-test analysis with Holm-Sidak correction. Retrosplenial: R; V1: Visual, Higher order visual: V+; Auditory: A; Somatosensory Secondary: Ss; Somatosensory Core: Sc; Somatosensory Barrel: Sb; Somatosensory other: S; Motor medial: Mm; Motor lateral: Ml.

Fig 10. Validation of domain maps by functional activity during locomotion.

(A) Behavioral video was used to extract locomotion/treadmill (red) and whisking (orange) time course. Boxed regions show analysis via optic flow previously described [41] (B) Correlation domain maps shown with either the re-addition of the original mean (left) and the filtered mean (right). Scatter plot showing the fraction of time walking or whisking (left) and correlation of the global mean to locomotion (right). (C) Comparison of the original mean (top; blue) and filtered mean (black) to the treadmill (red) and whisking (orange) time courses, and individual domain time courses with either the filtered (black) or original (blue) mean re-added that correlated to locomotion or were anticorrelated/less correlated. (D) Average domain map with global mean (left) and filtered mean (right) readded, functionally discrete regions highlighted. (E) Cluster analysis on the full movie (left) or just the times when locomotion was occurring (right). The co-clustering of the frontal motor region with barrel fields occurred in 6/7 mice.