2.1.

Animals

Four 3- to 4-month-old mice (two male, two female) were imaged at baseline (Day 0) and on Day 3 post photothrombotic stroke to left somatosensory forepaw cortex. All mice were Thy1-GCaMP6f [Jackson Laboratories Strain: C57BL/6J-Tg(Thy1-GCaMP6f)GP5.5Dkim; stock: 024276]. These mice express the protein GCaMP6f in excitatory neurons, primarily in cortical layers ii, iii, v, and vi.⁸ All studies were approved by the Washington University School of Medicine Animals Studies Committee and follow the guidelines of the National Institutes of Health’s Guide for the Care and Use of Laboratory Animals.

2.2.

Surgical Preparations

Prior to imaging, typical surgical preparations were implemented.^9,17 Briefly, isoflurane anesthesia (3% induction, 1% maintenance, $0.5\mathrm{L}/\min$) was used for sedation and an optically transparent $14\times 18\mathrm{mm}$ plexiglass window was implanted with translucent dental cement (C&B-Metabond, Parkell Inc., Edgewood, New York) following a midline incision and clearing of skin and periosteal membranes. The window covered the majority of the dorsal cortical surface and provided an anchor for head fixation and allowed for chronic, repeatable imaging.

2.3.

Photothrombosis

Mice were secured in a stereotaxic frame under isoflurane anesthesia. $200\mu \mathrm{L}$ of Rose Bengal (Sigma Aldrich) dissolved in saline ($10\mathrm{g}/\mathrm{L}$) was injected intraperitoneally. After 4 min, a 532-nm diode-pumped solid-state laser (Shanghai Laser & Optics Century) was focused to 2.2 mm left and 0.5 mm anterior to bregma with a 0.5 mm spot size and at 23 mW for 10 min.¹⁸ Mice were imaged at baseline [i.e., prior to photothrombosis (Day 0)], and 72 h post (Day 3). The dataset used in the following analyses consists of two five-minute imaging runs from each mouse. The stroke data were processed and analyzed as described below. Calcium data were filtered with a 0.4 to 4.0 Hz Butterworth bandpass filter and hemoglobin data with a 0.009 to 0.08 Hz Butterworth bandpass filter. These frequency bands were selected as they correspond to delta (0.4 to 4.0 Hz) and infraslow (0.009 to 0.08 Hz) ranges. The canonical FC frequency band (infraslow, 0.009 to 0.08 Hz) was used for hemoglobin-based analysis similar to the blood oxygen level dependent (BOLD) analysis used in the fMRI community.

2.4.

Fluorescence and Optical Intrinsic Signal (OIS) Imaging

Mice were head-fixed in a stereotaxic frame and body secured in a black felt pouch for imaging. Sequentially firing LEDs (Mightex Systems, Pleasanton, California) passed through a series of dichroic lenses (Semrock, Rochester, New York) into a liquid light guide (Mightex Systems, Pleasanton, California) that terminated in a 75 mm $f/1.8$ lens (Navitar, Rochester, New York) to focus the light onto the dorsal cortical surface. LEDs consisted of 470 nm (GCaMP6f excitation), 530, 590, and 625 nm light. An sCMOS camera (Zyla 5.5, Andor Technologies, Belfast, Northern Ireland, United Kingdom) coupled with an 85 mm $f/1.4$ camera lens (Rokinon, New York, New York) was used to capture fluorescence/reflectance produced at 16.8 Hz per wavelength of LED. A 515 nm longpass filter (Semrock, Rochester, New York) was used to discard GCaMP6f excitation light. Cross polarization (Adorama, New York, New York) between the illumination lens and collection lens discarded artifacts due to specular reflection. The field-of-view (FOV) recorded covered the majority of the convexity of the cerebral cortex ($\sim 1.1{\mathrm{cm}}^2$), extending from the olfactory bulb to the superior colliculus. All imaging data were binned in $156\times 156\text{pixel}$² images at $\sim 100{\mu \mathrm{m}}^2$ per pixel.

2.5.

Toolbox Capabilities and Workflow: Data and User Input

In order to initiate use of the toolbox, data has to be loaded into MATLAB (we recommend version 2022a or newer). Data should be in the form of pixels by pixels by frames. Example data used in the following analyses were acquired using the aforementioned mesoscopic calcium imaging modality, however, usage can be expanded to incorporate any data configured into this data stack (e.g., voltage-sensitive dye imaging, data from other animal models). However, one caveat to be noted is the following data processing pipeline was optimized on mesoscopic WOI mouse data. All inputs and outputs to specific scripts mentioned and highlighted in Fig. 1 are specified in the header of each script. The explanation that follows will walkthrough each processing step, which are all in separate subroutines. Subroutines will be referenced, and scripts that should be edited to run by the user are highlighted in Fig. 1 and Table 1.

Fig. 1

WOI mouse toolbox flowchart. (a) User input is needed after loading data, MATLAB will prompt user to trace the dorsal cortical surface as well as select anatomical landmarks. (b) Optical system dependent processing steps. These processing steps depend on the amount of ambient light present, natural fluctuations, and LEDs used for spectroscopy. (c) Optical system independent steps. (d) Optional steps such as affine transform if comparing multiple mice, temporal filtering, or quality control outputs. (e) Seed-wise FC analysis steps. (f) Bilateral FC analysis steps. (g) Stimulus activation analysis steps. (h) Cluster-based statistics steps used to find cluster size threshold in number of pixels.

Table 1

Data processing functions to be run by user. Functions outlined at each processing step in Figs. 1(a)–1(d).

Within the “START” folder at https://github.com/brierl/Mouse_WOI/tree/main/START, the folder “Proc” contains the first script that should be used to load data in [following the flowchart in Fig. 1(a)]. The script load_data.m allows the user to load the image stack and normalize one frame (variable “frame5”) that will be used in the next script make_mask_landmarks.m. Within make_mask_landmarks.m, the built-in MATLAB function roipoly.m is used to prompt the user to create a binary mask representing brain regions (i.e., 1) or nonbrain regions (i.e., 0) saved to the variable “isbrain.” The user is to click along the perimeter of the brain displayed in the frame and double click upon closing the loop to create the binary file [example in Fig. 2(a)]. This file will be used later and multiplied through the image stack so only pixels corresponding to brain is analyzed within the FOV. Next, the frame is passed to the subroutine MakeSeedsMouseSpace.m within make_mask_landmarks.m where the user is prompted to click on the anterior midline suture landmark and then the lambda landmark. From these two coordinates, a prototype of seed regions used for FC analysis [explained below and example in Fig. 2(b)] is displayed on the dorsal cortical surface (corresponding to motor, somatosensory, visual cortices etc.). A box will pop up prompting the user to specify if they are content with the seed prototype. A “yes” is appropriate if the seeds appear symmetrical across midline and extend across the entire dorsal cortical surface. A “no” will prompt the user to identify the aforementioned landmarks again with the hope of creating a more symmetrical seed prototype. These landmarks are later used to affine transform images to a common image space. Herein end the steps that rely on prompted user input [outlined in Fig. 1(a)].

Fig. 2

Example processing and analysis step outputs. (a) A traced dorsal cortical surface. (b) Seeds that are automatically plotted after selecting the anterior suture and lambda (marked by $X$). (c) Plots of pixels by time, the squared temporal derivative, or the RMS of the temporal derivative to visualize fluctuations in GCaMP6 signal during the experiment. (d) Example GCaMP6 seed-wise FC maps in one mouse, one run. (e) Example bilateral FC maps in one mouse, one run. (f) An overlaid cortical parcellation to visualize cortical functional regions. (g) A GCaMP6 activation trace using the same experimental paradigm as in Fig. 5(a) in one mouse and one experimental run. Signal was averaged across pixels in the activation area and vertical lines indicate the stimulus presence.

2.6.

Toolbox Capabilities and Workflow: Data Processing

The data processing functionality of the toolbox represents our standard image processing pipeline, which has incorporated methodologies across multiple WOI groups and is outlined continuing in Fig. 1(b). The first processing script that should be run is proc1_sys_dep.m and outlines the processing steps that are optical system dependent. First, a frame of ambient baseline light levels collected off of the imaging system (Dark.tif, will need to be replaced on each imaging system although ours is available on GitHub) is subtracted from the time series data in the subroutine subtract_dark_ois.m. Dark.tif was created by averaging dark frames collected over 1 min. Data is then spatially and temporally detrended using the built-in MATLAB detrend.m function within the subroutine detrend_ois.m to account for global artifactual fluctuations in the data by space and time. Examples of global fluctuations affecting the time-series light levels include photobleaching, LED current drift, and nonuniformities introduced by skull anatomy. Next, the processing stream diverges based on whether data is being processed into fluctuations in hemoglobin or fluorescence data. First, for hemoglobin data, channels that will be solved for fluctuations are sent to the subroutine procPixel.m where the logmean.m function takes the pixel-wise log-ratio of each frame (raw intensity data, $F_{\text{trace}}$, over mean intensity data),

The next data processing script is Proc2.m and outlines the steps that are optical system independent [Fig. 1(c)]. The first subroutine, smoothimage.m, spatially smooths the data with a $5\times 5$ Gaussian filter (1.3 pixel standard deviation). The global signal, computed by taking an average of all time traces within the FOV, gets regressed from the time series data to reduce global sources of variance such as motion and pulse, thus enhancing network correlation specificity (subroutine gsr.m).

Figure 1(d) outlines the optional processing steps. Affine.m allows for cross-mouse averaging by affine transforming the data to common Paxinos atlas space using the previously determined landmarks. The last processing script (OIS_GCaMP_Filter.m) allows for optional temporal filtering. Here, calcium data were filtered with a 0.4 to 4.0 Hz Butterworth bandpass filter and hemoglobin data with a 0.009 to 0.08 Hz Butterworth bandpass filter. These frequency bands were selected as they correspond to delta (0.4 to 4.0 Hz) and infraslow (0.009 to 0.08 Hz) ranges. The final processing script allows for an output of data quality (check_mvmt.m). This script outputs a pixels by time snapshot of contrast fluctuation [an example of this script run on GCaMP6 data in one run of one mouse in Fig. 2(c)]. The final panel is a calculation of the global variance in the temporal derivative, a surrogate measure of artifact tone as determined in human diffuse optical tomography.²⁰

2.7.

Toolbox Capabilities and Workflow: Analyses

Multiple types of analyses are supported by this toolbox and we walkthrough the two we use to analyze the stroke model in Fig. 1(e) and Table 2 (seed-based FC) and Fig. 1(f) and Table 3 (bilateral FC) post data processing. To initiate seed-based FC, folder “FC” has calc_fc.m, which takes in the processed data (e.g., “all_contrasts”), mask “isbrain,” and a seed set “seedcenter” to extract time traces within each seed region (subroutine P2strace.m) and calculate Pearson r correlation coefficients between each time trace and the remaining pixels in the FOV (subroutine strace2R.m). The next script, visualize_fc.m plots each seed-based FC map [example output in Fig. 2(d)].

Table 2

Seed-wise FC analysis functions to be run by user. Table of functions to calculate seed-wise FC maps as outlined in Fig. 1(e).

Table 3

Bilateral FC analysis functions to be run by user. Table of functions to calculate bilateral FC maps as outlined in Fig. 1(f).

For bilateral FC calculations, start with calc_bilateral.m within folder “BilatFC,” which requires the same data and mask inputs as calc_fc.m (note: this calculation does not need a seed set). The subroutine CalcRasterSeedsUsed.m outputs a matrix of left hand and right hand side pixel pairs that a Pearson r correlation value will be calculated between. Subroutines fcManySeed.m and CalcBilateral.m calculates the Pearson r value and organizes the r matrix into an image across the dorsal cortical surface. The next script, visualize_bilat.m, outputs a figure as shown in Fig. 2(e) of the bilateral FC maps for each contrast analyzed. A parcellated dorsal cortical surface is provided in Fig. 2(f) to orient users to the cortical regions within the FOV.

Additional analysis supported by the toolbox includes stim activation plotting and tracing [outlined in Fig. 1(g) and Table 4] in folder “Stims.” The user should run calc_stims.m with data (“all_contrasts”) and mask (“isbrain”) inputs. This script has commented variables to define the length of stimulus blocks so that the data can be rearranged into pixels by pixels by stimulus block length by number of blocks. Then, a mean frame (calculated by averaging together frames when the stimulus is off) is subtracted from the time-series data. The script plots an average frame from when the stimulus is on and thresholds the image at 80% max activation. The pixels surviving this threshold are averaged spatially, and the trace across averaged stimulus blocks is plotted [example trace for one mouse in Fig. 2(g)].

Table 4

Stimulus activation analysis function to be run by user. Function to calculate stimulus activation maps and traces as outlined in Fig. 1(g).

Outlined in Fig. 1(h) and Table 5 are the scripts used to perform the cluster size-based statistical thresholding (in folder “Stats”) described below. The user should run cluster_threshold.m, which requires data (“all_contrasts”), a mask (“isbrain”), and an alpha level [the overall family-wise error (FWE) rate e.g., $p=0.05$] as inputs. The script calls subroutine FWHM_ParDer.m, which estimates the FWHM based on a covariance matrix of the spatial partial derivatives of the data. Then, the threshold for cluster size significance is solved for and output.

Table 5

Cluster statistical analysis function to be run by user. Function to calculate cluster size for statistical significance as outlined in Fig. 1(h).

To enhance visualization, all of these maps can be plotted onto our cortical parcellation map (as in Fig. S2 in the Supplementary Material) as determined by White et al.,⁷ using the Parcellations.mat file (variable “Parcels”) located in the “START” folder.

2.8.

Electrical Hindpaw Stimulus Imaging

A stimulus driven protocol was used to further show the utility of the statistical correction method described. Two GCaMP6 mice were surgical prepared as previously described. Mice were anesthetized and microvascular clips were applied to the left hindpaw. An isolated pulse stimulator was used to deliver 2 Hz, $300\mu \mathrm{s}$, 0.5 mA electrical stimuli. Each imaging block consisted of 5 s rest, 5 s stimuli, and 10 s rest. Each mouse underwent 14 imaging blocks. GCaMP6 data were processed as described above and spatially downsampled to 78 pixels² and temporally downsampled to a framerate of 8.4 Hz. No temporal filtering was used and the pixel-wise average time trace was subtracted. Average activation maps were calculated by averaging across imaging blocks then across mice. A one-sample $t$-test was used to compare activation areas to the null hypothesis of no activation areas.

2.9.

Cluster-Based Thresholding

A challenge with analyzing the statistical significance in functional imaging is managing the multiple comparisons problem. Here, we used a cluster size-based method that leverages the spatial connection between pixels and credits large clusters as having more statistical significance than small clusters with the same peak $t$-value. More specifically, we used a cluster size-based thresholding method (determined by the cluster-size limit, $k_{\alpha }$) to analyze two-dimensional FC maps and to ensure the FWE rate did not exceed 5%. A random field theory (RFT) approach was adapted from the fMRI and DOT literature^15,16 to fit our two-dimensional data. Using RFT, we are able to approximate the expected ($E[]$) number of clusters ($m$) in an image at a given $z$-score threshold [$Z_t$, Fig. 3(a)]:

Fig. 3

Random field theory identifies cluster size threshold for pixel-wise $t$-maps on GCaMP6 bilateral FC data. (a) The relationship between the pixel-wise false positive rate (determined by $Z_t$, vertical dashed black line is $Z_t=3.09$, which corresponds to $p=0.001$) and the expected number of clusters to survive thresholding due to chance ($E[m]$). (b) The relationship between the pixel-wise false positive rate (determined by $Z_t$, horizontal dashed black line is $Z_{t }=3.09$, which corresponds to $p=0.001$), the family-wise error rate (here, 0.05) and cluster size needed for significance. (c) Top row: average ($N=4$) bilateral FC maps pre (left) and 3 days post (right) photothrombotic stroke to left somatosensory forepaw cortex (marked by hollow black circle). Middle row: pixel-wise paired $t$-test (left) and thresholded image for pixels with $p<0.05$ (right). Bottom row: thresholded image with FWE = 0.05 using the cluster size-based method (left) and image with Bonferroni correction for multiple comparisons (right).

We are able to then approximate the expected number of pixels ($N$) above $Z_t$ as

The probability that $n$ will exceed any threshold, $x$, can be modeled by the exponential function

Providing an equation (with all previously solved for values) that varies with the inverse of the squared FWHM. We only want clusters to survive thresholding with a family-wise error rate of 0.05 ($\alpha$), therefore, we want to find the pixel count threshold, $k_{\alpha }$, that would result in at least one cluster surviving threshold when there is no true significant difference, 5% of the time. Essentially, we are asking for a $k_{\alpha }$ where the following is true:

Assuming no clusters are a true positive. Which, mathematically speaking, is the same as solving for 1 minus the probability that no clusters have a pixel count above threshold $k_{\alpha }$ [Fig. 3(b)]:

After substituting Eq. (7) into Eq. (12) and solving for $k_{\alpha }$.

Pixel-wise $t$-tests were performed to compare Day 0 and Day 3 stroke FC maps, and to compare stimulus activated maps to the null hypothesis of no activation. These $t$-test maps were then thresholded, leaving pixels with a $t$-value corresponding to $p<0.05$. If the number of pixels in the remaining clusters exceeded $k_{\alpha }$, these clusters were considered to have a statistically significant finding within them by the cluster size-based technique.

2.10.

Data and Code Availability

To promote validation and comparative analyses by external groups, data and specific code will be made available through requests to the corresponding author. The toolbox presented here, (for WOI processing and analysis) is available online ( https://github.com/brierl/Mouse_WOI/tree/main/START/) and was used for all present analysis. Stroke and stim data is available on Mendeley Data (stroke²² and stim²³) and detailed in Table 6.

Table 6

Sample data available on Mendeley Data.22,23 Rows 1–4: table of the sample unprocessed data available online corresponding to two, 5-min baseline and two, 5-min post stroke runs. Rows 5 and 6: sample processed 5-min stim data runs in 14 blocks of 20 s each broken down into 5 s of rest, 5 s of 2 Hz stims, and 10 s of rest.

3.1.

Toolbox Overview

We seek to distribute a comprehensive, easy-to-use, open-source toolbox for mouse WOI data processing and analysis. Within this toolbox, we have translated multiple techniques from the human fMRI literature to the WOI mouse world. The toolbox is largely split up into two pipelines (Fig. 1, Tables 1–5), data processing and postprocessing analysis. A walkthrough is provided in Sec. 2, but briefly, prompted user input is needed while running scripts to load data and mark brain regions and anatomical landmarks [Figs. 1(a), 2(a), and 2(b)]. The remaining processing steps outlined in Figs. 1(b)–1(d) are split up into optical system dependent (for spectroscopy and fluorescence data processing) and independent steps, as well as optional steps [Fig. 2(c)]. After processing the data, multiple analysis pipelines are supported. Figures 1(e)–1(h) outlines the steps for the analyses shown in later figures, namely FC calculations [Figs. 2(d) and 2(e)], stimulus activation plotting [Figs. 2(f) and 2(g)], and statistical testing.

3.2.

Cluster Size-Based Thresholding Applied to Stroke Data

The cluster size-based statistical thresholding method^15,16 is detailed in the Methods and was adapted from the fMRI and DOT literature. The method is used here to select clusters of size $k_{\alpha }$ (expressed in number of pixels) with $p<0.05$ by the pixel-wise paired $t$-test method that satisfies the pixel-wise false positive rate (set by $Z_t$) and overall family-wise error rate [Figs. 3(a) and 3(b)]. Using this cluster size cutoff, we were able to compare bilateral FC maps at baseline ($N=4$, Day 0) and 72 h post ($N=4$, Day 3) photothrombotic stroke [Fig. 3(c)]. Photothrombosis was induced in left somatosensory cortex, which resulted in loss of homotopic FC at Day 3. A pixel-wise $t$-test was performed and thresholded to only display regions with $p<0.05$ (note, this map is not corrected for multiple comparisons). Using the cluster-size based threshold (FWE = 0.05), we were able to localize a somatosensory anchored deficit. This deficit was overlaid and contained completely within the parcellated somatosensory cortex (Fig. S2 in the Supplementary Material). Due to the symmetrical nature of the bilateral FC calculation, loss of homotopic FC is visualized by a symmetrical deficit [green area in Fig. 3(c) upper right panel]. Using the Bonferroni correction for multiple comparisons, no regions survived this stringent cutoff, resulting in no significant differences between day 0 and day 3 with this method.

Further exploring cross-network alterations in the context of this stroke model, we calculated seed-wise FC maps at Day 0 and Day 3 (Fig. 4). Using the same statistical framework, we compared the results at Day 0 and Day 3 and illustrated significant network alterations bilaterally in the seed-wise somatosensory (Ss.) map and an unaffected visual cortex outside of the area of photothrombosis.

Fig. 4

Cluster size thresholding yields Ss based deficit and unaffected visual cortex using a seed-wise FC methodology with calcium data. Average ($N=4$) seed-wise FC maps using seeds in somatosensory (Ss.) and visual (Vis.) cortices at (first column) baseline and (second column) 3 days post photothrombotic stroke to left somatosensory forepaw cortex (marked by hollow black circle). Third column: pixel-wise paired $t$-test between Day 0 and Day 3. Fourth column: thresholded images for pixels with $p<0.05$. Fifth column: thresholded images showing regions of statistically significant change using a clustered-based methodology with FWE = 0.05. Sixth column: no pixels survived the Bonferroni correction for multiple comparisons.

Repeating the same analysis and using oxygenated hemoglobin as a contrast, we see a similar deficit within somatosensory cortex in the seed-wise FC maps (Fig. S1 in the Supplementary Material). However, with hemoglobin as a contrast, we lose statistical significance in this experiment as the somatosensory-based deficit does not survive cluster-based statistical thresholding. Additionally, regions of cortex survive our cluster-based statistical analysis when using a visual seed contrary to when this analysis was performed using calcium.

3.3.

Cluster Size-based Thresholding Applied to Stimulus Data

We imaged two GCaMP6 mice that underwent anesthetized electrical hindpaw stimuli experiments. The protocol consisted of multiple blocks of 5 s of rest, 5 s of 2 Hz stimulus, and 10 s of rest [Fig. 5(a)]. The average activation while the stimulus was delivered was mapped [Fig. 5(b)]. A one-sample $t$-test was thresholded for $p<0.05$ leaving multiple areas of cortex. A Bonferroni correction for multiple comparisons removed all activation areas while the cluster-based correction method localized the somatosensory-based activation while removing spurious artifacts.

Fig. 5

Cluster size thresholding plots the GCaMP6 cortical activation in response to electrical sensory stimulus. (a) The experimental paradigm consisting of 5 s of rest, 5 s of stimulus on, 10 s of rest, that is repeated throughout the imaging duration. (b) Temporally averaged frame across durations of stimulus on (left) followed by a one-sample $t$-test thresholded for areas of $p<0.05$ and respective cluster-based and Bonferroni corrections.