Data-driven segmentation of cortical calcium dynamics

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.
Methods
Mice 

All animal studies were conducted in accordance with the University of California, Santa Cruz Office of Animal Research Oversight and Institutional Animal Care and Use Committee protocols. P21-22 Snap25 GCaMP6s transgenic mice (JAX: 025111), Cx3cr1 GFP (JAX: 005582), and Aldh1 GFP (MGI: 3843271) were maintained on a C57/Bl6 background in UCSCs mouse facilities. To identify Snap25 GCaMP expressing mice, a single common forward primer (5’-CCC AGT TGA GAT TGG AAA GTG-3’) was used in conjunction with either transgene specific reverse primer (5’-ACT TCG CAC AGG ATC CAA GA-3’; 230 band size) or control reverse primer (5’-CTG GTT TTG TTG GAA TCA GC-3’; 498 band size). The expression of this transgene resulted in pan-neuronal expression of GCaMP6s throughout the nervous system. To identify GFP expressing mice a forward (5’-CCT ACG GCG TGC AGT GCT TCA GC-3’) and reverse (5’-CGG CGA GCT GCA CGC TGC GTC CTC-3’; 400 band size) PCR amplification was used to identify which animals had the GFP transgene. At the end of each recording session, the animal was either euthanized or perfused and the brain dissected.

Surgical procedure

All mice were anesthetized with isoflurane (2.5% in pure oxygen) for the procedure. Body temperature was maintained at 30C for the duration of the surgery and recovery using a feedback-regulated heating pad. Lidocaine (1%) was applied subcutaneously in the scalp, followed by careful removal of skin above the skull. Ophthalmic ointment was used to protect the eyes during the surgery. The cranium was attached to two head bars using cyanoacrylate, one across the occipital bone of the skull and the other on the lateral parietal bone. After the surgery was complete, mice were transferred to a rotating disk for the duration of the recording which started 1 hour after the last isoflurane was administered. At the end of the recording session, the animal was euthanized and the brain dissected.

Recording calcium dynamics

In-vivo wide-field fluorescence recordings were collected in a minimally invasive manner. Imaging through the skull by single-photon excitation light from two blue LED lights (470 nm; Thorlabs M470L3) produces a green fluorescent signal that is collected through coupled 50mm Nikon lenses (f5.6 / f1.2, optical magnification ∼ 1x) into a scientific CMOS camera (PCO Edge 5.5MP, 6.5µm pixel resolution). The top lens in the tandem lens system was used to focus on the cortical surface, thereby lowering the magnification slightly; anatomical representation for each pixel corresponded to 6.9±0.2µm (min: 6.7µm, max: 7.2µm). Excitation light was filtered with a 480/30 nm bandpass (Chroma Technology AT480/30x) and the emission signal was filtered with 520/36 nm bandpass (Edmund Optics 67-044). Data collection was performed in a dark, quiet room with minimal changes in ambient light or sound, thus the brain activity recorded was at resting state without direct stimulation. Raw data was written directly as a set of 16-bit multi-image TIFF files. Concurrent body camera recording was taken using infrared (IR) raspberry pi cameras (Raspberry Pi NoIR V2), recorded at 30 fps synced to the CMOS.

The total amount of data recorded for each animal was generally at least 40 min and the amount of time in between video segments was less than 1 minute. All analyzed data consisted of two sequential full spatial resolution recordings concatenated together giving 20 min of data sampled at 10 frames per second for each video decomposition.

All video segments consisted of a set of continuously collected images at 10 or 20 frames per second for 10 minutes. The total amount of recorded data for each animal was generally at least 40 min and the amount of time in between video segments was less than 1 minute. When more than 10 minutes of video data were used for single decompositions, multiple videos were concatenated together. All analyses were conducted on data recorded at 10Hz, except exploration of effects of resolution on data quality.

Spatial resolution analyses were performed on a single 10 minute recording at 10Hz. Spatial down sampling was conducted by taking the mean between groups of pixels in a dxd square, where d is the integer down sampling factor. Temporal resolution analyses were performed on a single 10 minute recording at 20Hz. This data was temporally down sampled by taking the mean between an integer number of subsequent frames.

ICA decomposition and saving

ICA was performed using FastICA, implemented through python’s sklearn decomposition. The ICA decomposition was applied to the spatially flattened (xy,t) 2-D representation of the video data under the cortical ROI mask. The mean time series is pre-subtracted from the array before SVD decomposition or ICA decomposition, since ICA cannot separate sources with a mean signal effect. The filtered, unfiltered mean, ICA components, mixing matrix, and associated metadata are all saved. Data is stored and saved in this flattened format for storage optimization. Components are locally spatially reconstructed for visualization in the GUI.

Requesting the full number of components resulted in extremely lengthy processing times. To reduce the processing time, the data was preprocessed through Singular Value Decomposition (SVD) whitening, and noise components were cropped. To ensure that no signal was lost, and there were ample dimensions left for ICA separation, the inflection point between SVD signal and noise floor was identified, and SVD components were reduced to include components equal to 5 times the SVD signal to noise cutoff value. This cutoff multiplier can be adjusted while ICA projecting.

After calculating and sorting the ICA results, excessive noise components are removed from the dataset for compression. The cutoff was determined by identifying the inflection point in the lag-1 autocorrelation distribution with a two-peaked KDE fit. Components were saved such that 75% of the components saved were signal or artifact, and 25% of the components saved were noise associated. If not enough noise components were returned by the ICA decomposition, there is a risk that signals were not sufficiently unmixed, so the ICA decomposition was repeated with a higher SVD cutoff until enough additional noise components were included. This was not necessary for any of our experiments in this paper, but can be useful in ensuring signal capture in recordings of smaller areas, such as cranial windows.

For resolution down sampling analysis, the number of components requested was kept constant for all resolution values to avoid confounding the results. The cutoff for the highest resolution tested was used for all lower resolution decompositions.

ICA returns components that are unsorted and often spatiotemporally inverted. Components were first sorted by their time series variance. The spatial histogram of sIC values across each component can be visualized as a single-tailed gaussian distribution centered around 0, where the tail represents the spatial domain of each component. The two edges of the distribution are first identified. The boundary closer to 0 is taken as the edge of the central noise distribution, and that boundary is used to define the dynamic threshold on the opposite side of 0. Any values outside of this noise distribution and is part of the wider tail are included in the binarized domain of the component. If the tail was negative, the component was inverted spatially and temporally by multiplying each by a factor of -1. In this way, components were all identified as positive effectors for visualization purposes. Movie rebuilding is not affected by this process.