Yujian Diao
Ting Yin
Rolf Gruetter
Ileana O. Jelescu- 1Animal Imaging and Technology, EPFL, Lausanne, Switzerland
- 2CIBM Center for Biomedical Imaging, Lausanne, Switzerland
- 3Laboratoire d’Imagerie Fonctionnelle et Métabolique, EPFL, Lausanne, Switzerland
Resting state functional MRI (rs-fMRI) is a widespread and powerful tool for investigating functional connectivity (FC) and brain disorders. However, FC analysis can be seriously affected by random and structured noise from non-neural sources, such as physiology. Thus, it is essential to first reduce thermal noise and then correctly identify and remove non-neural artifacts from rs-fMRI signals through optimized data processing methods. However, existing tools that correct for these effects have been developed for human brain and are not readily transposable to rat data. Therefore, the aim of the present study was to establish a data processing pipeline that can robustly remove random and structured noise from rat rs-fMRI data. It includes a novel denoising approach based on the Marchenko-Pastur Principal Component Analysis (MP-PCA) method, FMRIB’s ICA-based Xnoiseifier (FIX) for automatic artifact classification and cleaning, and global signal regression (GSR). Our results show that: (I) MP-PCA denoising substantially improves the temporal signal-to-noise ratio, (II) the pre-trained FIX classifier achieves a high accuracy in artifact classification, and (III) both independent component analysis (ICA) cleaning and GSR are essential steps in correcting for possible artifacts and minimizing the within-group variability in control animals while maintaining typical connectivity patterns. Reduced within-group variability also facilitates the exploration of potential between-group FC changes, as illustrated here in a rat model of sporadic Alzheimer’s disease.
Introduction
Resting state functional MRI (rs-fMRI) based on spontaneous low-frequency fluctuations in the blood oxygen level dependent (BOLD) signal in the resting brain is a widely used non-invasive tool for studying intrinsic functional organization in health and disease (Fox and Raichle, 2007; Fornito and Bullmore, 2010). By examining spatio-temporal correlations of the BOLD signal between distinct brain regions, known as functional connectivity (FC), this technique is capable of revealing large-scale resting state networks (RSNs) (Biswal et al., 1995; Damoiseaux et al., 2006; Buckner et al., 2013). Nowadays, rs-fMRI has become an increasingly important translational neuroimaging tool for understanding neurological and psychiatric diseases and for developing treatments, with rapidly growing applications not only in human research but also in rodent models of disease (Fox and Greicius, 2010; Bajic et al., 2017).
However, the BOLD signal is contaminated by multiple physiological and non-physiological sources of noise, such as respiratory and cardiac cycles, thermal noise, changes in blood pressure, and head motion (Kruger and Glover, 2001; Birn, 2012; Van Dijk et al., 2012; Murphy et al., 2013). These non-neuronal sources can severely affect rs-fMRI time series and thereby confound the connectivity analysis (Cole et al., 2010; Power et al., 2014). Therefore, a robust pre-processing pipeline is required to extract the neuronal component of the BOLD signal and minimize the contribution of such noise sources. Furthermore, existing tools that correct for the effect of non-neuronal sources are mostly tailored for human rs-fMRI data and are not readily transposable, or even applicable, to rodent data. Dedicated pipelines for rodent rs-fMRI processing are just starting to emerge (Zerbi et al., 2015; Bajic et al., 2017).
For example, signal fluctuations resulting from respiratory and cardiac cycles can be accounted for by explicitly recording the physiology and modeling these external confounds as regressors (Birn et al., 2006; Kasper et al., 2017). While physiological recordings in rodents are possible, they typically involve dedicated hardware and invasive procedures, making them experimentally difficult. However, although cardiac and respiratory frequencies in rodents are much higher than those of the resting-state BOLD fluctuations, depending on the temporal resolution of the acquisition, they can alias into the band of interest (typically 0.01–0.3 Hz) and corrupt the analysis. Two complementary approaches are, therefore, suitable to mitigate the impact of physiological noise in rodent rs-fMRI.
One approach is the removal of global signal defined as the mean time series averaged over all voxels within the brain by including the global signal as a nuisance regressor in General Linear Model (GLM) analyses, which is referred to as global signal regression (GSR) (Liu et al., 2017). However, the use of GSR has been one of the most controversial topics in human rs-fMRI connectivity studies (Liu et al., 2017; Murphy and Fox, 2017). On one hand, GSR is known to introduce spurious negative correlations (Murphy et al., 2009) and cause spatial bias on connectivity measures (Saad et al., 2012). On the other hand, prior studies have shown that GSR can enhance the detection of significant FC and improve spatial specificity of positive correlations (Fox et al., 2009). Most importantly for rodent studies, GSR can also mitigate confounds related to motion and physiological processes (Power et al., 2015; Aquino et al., 2019).
Another commonly used data-driven approach that identifies various physiological noise components directly from the fMRI data itself is single-level independent component analysis (ICA) (McKeown et al., 2003; Griffanti et al., 2014; Bajic et al., 2017; Caballero-Gaudes and Reynolds, 2017). The ICA method is also confronted by several issues including model order selection (i.e., the number of components) (Kuang et al., 2018) and the identification of artifactual components, which is a manually tedious step (Wang and Li, 2015), especially for a high order model. Notably, a machine-learning approach for automatic artifact component classification based on FMRIB’s ICA-based Xnoiseifier (FIX) (Salimi-Khorshidi et al., 2014) has been proposed to replace manual classification. The FIX auto-classifier applied in human and mice rs-fMRI studies has yielded promising results with a high accuracy in artifact identification (Griffanti et al., 2014, 2015; Zerbi et al., 2015). However, the success of FIX classification relies on a proper pre-training on study-specific datasets.
Therefore, the aim of the present study was to propose and evaluate a data processing pipeline for rat rs-fMRI that minimizes intra-group variability and maximizes between-group differences in whole-brain FC. In this pipeline, we reduced structural noise by combining single-session ICA cleaning and GSR. For ICA cleaning, we built and used a dedicated FIX classifier for rats. Furthermore, we enhanced the sensitivity to the BOLD fluctuations by first increasing dramatically the temporal signal-to-noise ratio (tSNR) of the data. For the purpose of stochastic (thermal) noise removal, we employed a novel method based on Marchenko-Pastur Principal Component Analysis (MP-PCA). MP-PCA denoising was recently introduced for diffusion MRI and is a model-free method that exploits redundancy in MRI series (Veraart et al., 2016), which has shown great potential for improving the SNR in other MRI techniques as well (Does et al., 2019; Ades-Aron et al., 2020).
Materials and Methods
Animal Preparation and Anesthesia
All experiments were approved by the local Service for Veterinary Affairs. Male Wistar rats (236 ± 11 g) underwent a bilateral intracerebroventricular (icv) injection of either streptozotocin (3 mg/kg, STZ group) or buffer (control group). When delivered exclusively to the brain, streptozotocin induces impaired brain glucose metabolism and is used as a model of sporadic Alzheimer’s disease (AD) (Lester-Coll et al., 2006; Knezovic et al., 2015; Grieb, 2016).
As a result of system upgrade, the fMRI data were acquired on two rat cohorts (1/2) on different MRI consoles (Varian/Bruker). We refer to animals scanned on the Varian system as Cohort 1 (N = 17 rats) and to animals scanned on the Bruker system as Cohort 2 (N = 7 rats). Each cohort comprised animals from both groups: Cohort 1 (CTL/STZ, N = 8/9 rats) and Cohort 2 (CTL/STZ, N = 4/3 rats), which were scanned at 2, 6, 13, and 21 weeks after the surgery (Figure 1). Rats were anesthetized using 2% isoflurane in a mixture of O2 and air (O2/air: 30/70) during the initial setup and promptly switched to medetomidine sedation delivered through a subcutaneous catheter in the back (bolus: 0.1 mg/kg, perfusion: 0.1 mg/kg/h) as previously described (Reynaud et al., 2019). Medetomidine preserves neural activity and vascular response better than isoflurane (Weber et al., 2006; Pawela et al., 2009; Grandjean et al., 2014). The rat head was fixed using a homemade holder with a bite bar and ear bars to minimize the head motion, and body temperature and breathing rate were continuously monitored. At the end of the scanning sessions, rats were woken up with an intra-muscular injection of antagonist atipamezole (0.5 mg/kg) and returned to their cages.
Figure 1. Timeline of experiments. Two fMRI runs were acquired per rat for each experiment. Figure adapted from Tristão Pereira et al. (2021), NeuroImage 2021, with permission.
MRI Acquisition
Magnetic resonance imaging experiments were conducted on a 14.1 T small animal scanner with two different consoles: Varian system (Varian Inc.) equipped with 400 mT/m gradients (Cohort 1) and Bruker system (Ettlingen, Germany) equipped with 1 T/m gradients (Cohort 2), both using an in-house built quadrature surface transceiver. The acquisition parameters were the same for the two cohorts.
An anatomical reference scan was acquired using a fast spin-echo multi-slice sequence with the following parameters: TE/TR = 10.17/3000 ms, ETL = 4, TEeff = 10.17 ms, field of view (FOV) = 19.2 mm × 19.2 mm, matrix = 128 × 128, in-plane resolution = 150 μm × 150 μm, number of slices = 30, and thickness = 0.5 mm. Before running the fMRI sequence, anesthesia was switched from isoflurane to medetomidine. The fMRI acquisitions were started after a fixed duration (∼1 h) since the switch from isoflurane to medetomidine to minimize between-animal anesthesia-related confounds. Rs-fMRI data were acquired using a two-shot gradient-echo echo-planar imaging (EPI) sequence as follows: TE = 10 ms, TR = 800 ms, TRvol = 1.6 s, FOV = 23 mm × 23 mm, matrix = 64 × 64, in-plane resolution = 360 μm × 360 μm, 8 slices, thickness = 1.12 mm, and 370 repetitions (∼10 min). Two fMRI runs were acquired for each rat. A short scan with 10 repetitions and reversed phase-encoding direction was also acquired to correct for EPI-related geometric distortions.
FMRI Data Pre-processing
Our data processing pipeline was developed based on the data of Cohort 1.
Anatomical and fMRI images were first skull-stripped separately using BET (Brain Extraction Tool; FSL)1 (Smith, 2002), and fMRI time series were denoised using MP-PCA with a 5 × 5 × 5 voxel sliding kernel (Veraart et al., 2016). The quality of MP-PCA denoising was assessed by inspecting the normality of the residuals (original – denoised) and the tSNR changes before and after denoising. Specifically, the normality of the residuals was tested by the linearity of the relationship between the natural logarithm of the residual distribution probability and the squares of multiple residual standard deviation. Then, the datasets went through EPI distortion correction using FSL’s topup (Smith et al., 2004), slice-timing correction (Henson et al., 1999; Calhoun et al., 2000; Sladky et al., 2011), and spatial smoothing (Gaussian kernel: 0.36 mm × 0.36 mm × 1 mm). Corrected fMRI images were registered to the Waxholm Space Atlas of the rat brain2 using linear and non-linear registration in ANTs (Avants et al., 2008). Highly parcelated atlas labels were grouped to yield larger, anatomically consistent, labels [e.g., all sub-regions pertaining to the anterior cingulate cortex (ACC) were grouped into one ACC label], and 28 atlas-defined regions of interest (ROIs, 14 per hemisphere) were thus automatically segmented. The brain mask extracted based on EPI images was used to mask atlas labels in fMRI space such that only regions with sufficient signal were included and areas of drop-out were excluded.
Finally, single-session ICA was performed on fMRI time courses using FSL’s MELODIC (Beckmann and Smith, 2004) with high-pass temporal filtering (f > 0.01 Hz) and 40 independent components (ICs).
FIX Training
Datasets of Cohort 1 were randomly split into two groups: a training dataset for FIX (n = 49) and a test dataset (n = 58). The ICA components in the training dataset were manually classified to signal or artifact, which was mainly based on thresholded spatial maps because ICA is theoretically more robust in the spatial than in the temporal domain (Smith et al., 2012; Salimi-Khorshidi et al., 2014). Generally speaking, spatial maps of signal components should contain a low number of anatomically consistent clusters, whereas artifactual components typically have either very large clusters covering brain slices or very small and scattered clusters (Griffanti et al., 2017). Here, we chose an “aggressive” artifact removal (Griffanti et al., 2014) in the training dataset in order to give the trained classifier a margin to be conservative or aggressive via adjusting the threshold fed to it (small thresholds make it conservative).
The performance of the trained classifier in detecting artifactual components was evaluated on the test dataset by comparing the automatic classification of artifactual components with the manual classification. The classification accuracy was characterized in terms of “recall” and “precision” (Powers, 2011), which are defined as the percentage of the correctly predicted artifact components in all actual artifact components and the percentage of the correctly predicted artifact components in all predicted artifact components, respectively.
Network Analysis and GSR
After ICA decomposition and classification, the artifactual components were regressed out of the 4D pre-processed datasets to obtain “cleaned” rs-fMRI datasets. The cleaned data were used to compute ROI-to-ROI FC by calculating correlation coefficients between the ROI-averaged time series of the 28 atlas-defined ROIs, resulting in a 28 × 28 FC matrix for each rat.
For FIX, we relied on spatial maps rather than time courses to evaluate artifactual components. In several cases, components displayed sensible spatial distribution to represent a RSN, but the power spectrum showed a peak at a frequency that could be attributed to breathing (Supplementary Figure 1) (the breathing rate was recorded for each run and its aliased frequency within our 0.01–0.31 Hz band was calculated). These non-neuronal sources that could not be cleaned using FIX also contribute to the global signal. To mitigate their effect, we used the partial correlation of ROI-to-ROI time courses to build FC matrices, with the global signal as the controlling variable. For every pair of ROIs, the partial correlation was implemented by measuring the correlation between their time-series residuals, after each having been adjusted by the GSR (Smith et al., 2011).
Finally, statistical comparisons of FC between the STZ and CTL groups at each timepoint were performed using NBS (Zalesky et al., 2010) to identify network connections that showed significant between-group difference. Specifically, NBS uses one-tailed two-sample t-test to detect differences in group-averaged FC between the two groups. Thereby, two contrasts (STZ > CTL and STZ < CTL) were tested separately. A t-statistic threshold was chosen on the basis of medium-to-large sizes of the subnetwork comprised connections with their t-statistic above the threshold (Tsurugizawa et al., 2019) as well as the underlying p-values. Here, we chose 2.2 as the t-statistic threshold. Results based on other thresholds are provided as Supplementary Figures 2–4. Significance (p ≤ 0.05) was tested after family wise error rate (FWER) correction using non-parametric permutation (N = 5000).
The full data processing pipeline is illustrated in Figure 2.
Figure 2. The proposed pipeline for rs-fMRI data processing.
Pipeline Evaluation on Cohort 1
To evaluate the rs-fMRI data processing pipeline including MP-PCA denoising (DN), slice-timing correction (SC), spatial smoothing (SM), high-pass filtering (HP), ICA-FIX cleaning (CL), and GSR in terms of consistency of within-group FC in the healthy CTL group and in terms of between-group difference, we compared results of the PIRACY (Pipeline for Rat Connectivity) pipeline (DN + CL + GSR, pipeline G) with other processing approaches depending on the presence of DN, CL, and GSR in addition to the baseline pipeline (SC + SM + HP, pipeline A), namely, GSR (pipeline B), CL + GSR (pipeline C), DN (pipeline D), DN + CL (pipeline E), and DN + GSR (pipeline F) in Cohort 1, shown in Table 1.
Table 1. The seven data processing pipelines and methods they include (“×” – including, “∘” – excluding).
Based on the hypothesis that an optimal processing procedure should minimize the variability within the homogeneous group of healthy controls (Zerbi et al., 2015), the within-group variability was assessed by calculating the standard deviation of the Fisher z-transformed correlation coefficients of the FC matrices in the CTL group of healthy rats at each timepoint. In addition, the sensitivity to between-group differences was evaluated by comparing the significant difference in FC between the STZ and CTL groups at each timepoint.
Group ICA Analysis
At each timepoint, group-level ICA was performed using FSL’s Melodic (Smith et al., 2004) on the pooled CTL and STZ datasets of Cohort 1 that were pre-processed using pipeline E (no GSR). Prior to group-level ICA in Melodic, all the rs-fMRI datasets from a given timepoint were registered to a common template using ANTs (Avants et al., 2008), and the template registration in Melodic was by-passed. Thirty group-level spatial ICs were extracted. The selection of IC number is a trade-off between biological detail (small and refined networks) and noise (too many components that can be supported by data quality) (Bajic et al., 2017). N = 30 or 40 components are common for rodent group ICA (Jonckers et al., 2011; Zerbi et al., 2015; Bajic et al., 2017). In our case, N = 30 had the advantage of yielding a similar number of nodes to the ROI-to-ROI FC analysis for the network analyses. Dual regression (Filippini et al., 2009) was used to estimate subject-specific time courses and associated spatial maps. Similar to the ROI-based FC analysis, ICA-based FC matrices were built by calculating correlation coefficients between all pairs of time courses (excluding artifact components) using FSLNets (Griffanti et al., 2014). Both full correlation (original Pearson’s correlation) and partial correlation with GSR (mean time course as regressor) were employed by FSLNets to evaluate the similarity between time courses. In the end, differences in FC between the STZ and CTL groups were tested using the NBS toolbox using the same parameters in the Network Analysis and GSR section.
Pipeline Evaluation on Cohort 2
To test the robustness of the proposed pipeline, we further evaluated it on the independent datasets of Cohort 2. Similarly, we compared results of the PIRACY pipeline G (DN + CL + GSR) with three other processing pipelines: DN (D), DN + CL (E), and DN + GSR (F) in terms of the within-group variability for the control rats of Cohort 2. Then, we assessed the between-group differences in FC between the STZ and CTL groups on the pooled datasets of Cohorts 1 and 2.
Results
To experimentally evaluate the processing pipeline, a total of 109 rs-fMRI datasets (Cohort 1) were acquired from 17 rats at four timepoints ranging from 2 to 21 weeks (Figure 1). Two datasets with bad image quality were discarded. Figure 3 shows the MR images in one representative dataset including rs-fMRI images, matching anatomical reference, and the atlas-based anatomical labels registered to the fMRI images.
Figure 3. Example of rs-fMRI images of eight coronal slices (A), matching anatomical MR images (B), and atlas-based anatomical labels registered to the fMRI images (C).
To assess the robustness of the PIRACY pipeline, an independent cohort with 56 rs-fMRI datasets was acquired from seven different rats at four timepoints ranging from 2 to 21 weeks (Figure 1).
MP-PCA Denoising
The average tSNR after MP-PCA denoising improved significantly for all the 107 datasets. Figure 4 shows an example of the average tSNR increase from 75 to 146 after MP-PCA denoising, residuals map, histogram, and the normality test. The linearity of log(P) = f(r2) confirms that the residuals are normally distributed and only Gaussian noise has been removed from the signal. Moreover, the number of components classified as artifacts decreased, whereas the z-statistic for signal components increased when MP-PCA denoising was applied prior to ICA decomposition (Figure 5).
Figure 4. Example of temporal SNR maps, before (A) and after MP-PCA denoising (B). The mean tSNR over the middle brain slice was improved dramatically from 75 to 146. The SNR profile is typical of a surface coil placed on top of the head, with higher sensitivity in the cortex. (C) Residuals map of all voxels within the brain mask (rows) and time frame (columns). (D) Histogram of residuals. (E) Normality test.
Figure 5. (A) Number of ICs classified as artifacts across all datasets without (blue) and with (red) prior MP-PCA denoising, respectively. MP-PCA denoising substantially reduced artifactual ICA components for most datasets. (B) Prior MP-PCA denoising improved the z-statistic for signal components. Thresholded IC maps in (B) were given by FSL’s MELODIC.
FIX Classification
Here, we preferred a lower order model to avoid overfitting (Li et al., 2007; Kuang et al., 2018), and we chose the number of ICs to be 40, which typically explained 90% of the variance. Notably, reaching 95% of explained variance would have required about 90 components, which would potentially cause over-splitting networks and making the classification more complicated.
The single-subject ICA was performed on each dataset with 40 components. In the training dataset (n = 49), 19.8 ± 4.7 components (50%) were classified as artifacts by hand. In the test dataset (n = 58), between 40 and 64% of components were classified as artifacts automatically by FIX depending on the threshold. More components could be recognized as artifacts by increasing the FIX threshold at the expense of lowering the classification precision due to more misclassification. Here, 45 might be an “optimal” threshold with overall 88% in recall and 90% in precision achieved (Table 2). In practice, we ran FIX twice with a low and high threshold, respectively. We then examined only the components differing between the two artifact lists given by the two FIX thresholds and manually restored the true noise components. This represented nonetheless a significant gain in processing time, as we only needed to classify one third of the ICA components by hand.
Table 2. FIX artifact classification accuracies at different thresholds.
Comparison of Pipeline Performance
The seven data processing pipelines (A–G) were first compared based on resulting FC matrices of 28 atlas-based ROIs in Cohort 1.
Pipelines C and G that include both CL and GSR yielded lower within-group variability in the homogeneous group of healthy controls for all timepoints than other procedures excluding CL and/or GSR (Figure 6). The combination of CL + GSR also corrected for physiological artifacts that CL or GSR alone could not systematically address (Figure 7). Furthermore, our proposed pipeline G (DN + CL + GSR) yielded between-group differences most consistently at 2, 6, and 13 weeks, whereas pipeline C (without DN) only detected between-group differences at 2 and 6 weeks. No significant difference was found between the CTL and STZ groups for the baseline pipeline A and pipeline D (only DN). Pipeline E (DN + CL) exhibited numerous different edges at 6 weeks and none at the other timepoints, whereas pipelines B and F (GSR or DN + GSR) yielded numerous different edges at 6 and 21 weeks but none at 2 and 13 weeks (Figure 8).
Figure 6. Standard deviation of the Fisher z-transformed correlation coefficient for functional connectivity in the CTL group of Cohort 1 at four different timepoints for the seven pipelines (A–G). Pipelines including CL and GSR (C,G) obtained the minimal within-group variability in the homogeneous group of healthy controls for all timepoints, whereas other procedures excluding CL and/or GSR had higher variability.
Figure 7. Example FC matrices of two control rats scanned at 2 weeks (Subjects 1 and 2) processed with four different pipelines (D–G). For Subject 1, GSR alone (pipeline F) reinforces a widespread anti-correlation between the left and right hemispheres, manifest in a checkerboard pattern (odd/even rows and columns correspond to the left/right hemisphere ROIs, respectively) that was otherwise removed by ICA cleaning (pipeline E). For Subject 2, GSR alone does not attenuate very strong positive/negative correlation blocks spanning large portions of the brain, whereas ICA cleaning alone yields artificially high whole-brain connectivity. The combination of CL + GSR (pipeline G) mitigates these effects in both examples and yields more consistent between-subject connectivity matrices.
Figure 8. Significant differences in FC between the CTL and STZ groups of Cohort 1 at each timepoint (2, 6, 13, and 21 weeks) for each data processing approach under NBS threshold 2.2. No significant difference was found between the CTL and STZ groups for the baseline pipeline (A) and pipeline (D) (DN only), whereas significant differences were only found at 6 weeks for pipeline (E) (DN + CL, no GSR). Pipelines (B,F) (GSR, or DN + GSR) had similar results with between-group differences detected at 6 and 21 weeks. Pipeline (G) (DN + CL + GSR) most consistently yielded between-group differences across timepoints, whereas pipeline (C) (CL + GSR, no DN) only detected significant differences at two timepoints. Blue edges indicate group differences in contrast 1 (STZ > CTL), and red edges indicate group differences in contrast 2 (STZ < CTL). p1 and p2 are FWER corrected p-values for contrast 1 and contrast 2, respectively. Green nodes indicate ROIs involved in group differences.
ROI-Based FC Analysis
To analyze group differences in connectivity obtained following the PIRACY pipeline G, a more detailed visualization is provided in Figure 9. The complete list of edges with significant between-group differences is provided as Supplementary Table 1. Nodes with altered connectivity in STZ animals were consistent with regions affected by AD and with previous findings on this animal model (Mayer et al., 1990; Shoham et al., 2003; Lester-Coll et al., 2006; Kraska et al., 2012; Grieb, 2016). At 2 weeks, increased positive connectivity of the ACC to the retrosplenial cortex (RSC) and decreased anti-correlations of the default mode network (DMN) including the ACC, RSC, posterior parietal cortex (PPC), and hippocampus (Hip/Sub) to the lateral cortical network (LCN) involving somatosensory (S1) as well as motor (M) were found. The 6-week timepoint showed widespread reduced anti-correlations between the DMN (including the RSC, PPC, Hip), somatosensory of the LCN, and striatum (Str). At 13 weeks, weaker positive correlations were found within the DMN involving not only the ACC, PPC, RSC, Hip, medial temporal lobe (MTL), and visual cortex (V) but also the hypothalamus (HTh) to the DMN.
Figure 9. (A) Hybrid average FC matrices at each timepoint (top-right half: STZ, bottom-left half: CTL) based on the data of Cohort 1 processed by the optimized pipeline G. ∗: p < 0.05 (FWER corrected) at threshold of 2.2. (B) Graph networks at four timepoints. Blue/red edges and green nodes indicate connections with significant difference. The p-value for each network was given after FWER correction. ACC, anterior cingulate cortex; RSC, retrosplenial cortex; PPC, posterior parietal cortex; MTL, medial temporal lobe; Hip, hippocampus; Sub, subiculum; Au, auditory; V, visual; S1/S2, primary/secondary somatosensory; M, motor; Str, striatum; Tha, thalamus; HTh, hypothalamus; L/R, left/right. Figure adapted from Tristão Pereira et al. (2021), NeuroImage 2021, with permission.
Group ICA
Group ICA with 30 components was carried out on Cohort 1. Artifactual components were identified and removed, which resulted in 28, 25, 25, and 26 signal components left for the four timepoints, respectively. Figure 10 displays the significant between-group differences in partial correlations with GSR between ICA components. The complete list of edges with significant between-group differences is provided as Supplementary Table 2. Remarkably, group ICA analysis with GSR showed intergroup differences at three timepoints from 2 to 13 weeks, in agreement with differences found in ROI-based FC analysis using pipeline G. At 2 weeks, most changes were found in connections between the RSC, PPC, Hip, thalamus (Tha), and S1/2. At 6 weeks, there were alterations found in connectivity involving the PPC, Hip, RSC, S1/2, Str, as well as Tha. At 13 weeks, primary changes were detected in connections between the Hip, S1, Str, Tha, and HTh. Full correlation between ICA components–an equivalent to pipeline E: DN + CL–revealed significant differences at all timepoints (Supplementary Figure 5), in poor agreement with the ROI-based FC analysis using pipeline E, which only identified group differences at one timepoint (6 weeks). This inconsistency may indicate the importance of GSR in obtaining consistent intergroup differences.
Figure 10. Between-group differences in ICA-based functional connectivity with GSR (partial correlation) at NBS threshold of 2.2 for each timepoint in Cohort 1. Colored edges display the existence of STZ–CTL difference in connections between ICs. Each IC is denoted by a spatial map and its IC number. The nodes of ICs are listed in an order based on its position in the brain (anterior to posterior). ROI labels are attached to every ICA component. Artifactual components were removed, and the IC number was reordered accordingly. Blue edges indicate group differences in contrast 1 (STZ > CTL), and red edges indicate group differences in contrast 2 (STZ < CTL).
Pipeline Performance on Independent Datasets (Cohort 2)
The proposed pipeline G also yielded the minimal within-group variability in control rats of Cohort 2 at all timepoints, compared with three other pipeline variants without GSR and/or cleaning, in agreement with results on Cohort 1 (Figure 11). Moreover, between-group differences in FC in the combined datasetCohorts 1 + 2 (Figure 12) were consistent with that in Cohort 1 alone (Figure 6), in a pattern of initial hyperconnectivity and later hypoconnectivity from 2 to 13 weeks. At 21 weeks, Cohort 1 alone was likely underpowered due to missing datasets, and no group differences were reported with pipeline G. The combination of Cohorts 1 + 2 enabled to reveal widespread hypoconnectivity in STZ at 21 weeks (see Supplementary Figure 6 for a detailed comparison of significant edges for Cohort 1 and Cohorts 1 + 2 at all timepoints). In contrast, group differences as highlighted by pipelines E (DN + CL) and F (DN + GSR) showed less consistency between Cohort 1 and Cohorts 1 + 2, with dramatic changes in outcome at 6 weeks particularly.
Figure 11. Standard deviation of the Fisher z-transformed correlation coefficient for functional connectivity in the CTL group of Cohort 2 at four different timepoints for four pipelines (D–G) all including DN. Pipeline (G) (DN + CL + GSR) yielded the minimal within-group variability in the homogeneous group of healthy controls for all timepoints, whereas other procedures excluding CL and/or GSR had higher variability.
Figure 12. Significant differences in FC between the CTL and STZ groups of pooled Cohorts 1 + 2 at each timepoint (2, 6, 13, and 21 weeks) for data processing approaches (D–G) under NBS threshold 2.2. No significant difference was found with pipelines (D–F) (except pipeline F at 21 weeks). In particular, previous widespread differences at 6 weeks with pipelines (E,F) from Cohort 1 (Figure 8) disappeared entirely for Cohorts 1 + 2. Pipeline (G) exhibited similar group differences to previously (Figure 8) with the exception of new significant edges at 21 weeks, given a large increase in number of datasets at 21 weeks with the addition of Cohort 2. Blue edges indicate group differences in contrast 1 (STZ > CTL), and red edges indicate group differences in contrast 2 (STZ < CTL). p1 and p2 are FWER corrected p-values for contrast 1 and contrast 2, respectively. Green nodes indicate ROIs involved in group differences.
Discussion
In this work, we proposed a novel rs-fMRI processing pipeline adapted for rat data: PIRACY, which combines MP-PCA denoising, ICA cleaning, and GSR. We built a dedicated FIX ICA classifier for rat brain that showed a high accuracy in distinguishing artifactual ICA components from the rs-fMRI signal after training. We evaluated the performance of the proposed pipeline by comparing it to six other possible approaches that excluded denoising, artifact cleaning, and/or GSR. We further evaluated the pipeline on separate datasets from an independent cohort. We showed that these three steps were essential in minimizing the within-group variability in the healthy control group. Differences between a control group and a diseased group using the STZ animal model were found more consistently across timepoints with pipeline G and were also more consistent between different analysis approaches: ROI-to-ROI FC or group-level ICA.
We included a novel thermal noise reduction method based on MP-PCA applied to rs-fMRI data that substantially improved the tSNR, resulting in fewer artifactual components to be cleaned in the ICA decomposition, and improved the statistical strength and anatomical consistency of components classified as signal. While the MP-PCA denoising technique was initially implemented for diffusion MRI data (Veraart et al., 2016), it has recently also shown great value for task fMRI (Ades-Aron et al., 2020), and it is also applicable to rs-fMRI. Indeed, resting-state neural activity may appear as random temporal events, but it falls into a specific 0.01–0.3 Hz frequency range, which makes it discernible from white Gaussian thermal noise, with a flat power spectral density. This feature is apparent in the normality test (Figure 4E) for the residuals of MP-PCA denoising. Furthermore, denoising is performed within a sliding window of 5 × 5 × 5 voxels, within which the BOLD fluctuations are likely to be correlated (cluster), whereas the thermal noise is not. Thus, the thermal noise reduction using the MP-PCA approach is very unlikely to remove the genuine BOLD fluctuations and, on the contrary, will improve the sensitivity of the analysis for the latter, as is apparent from the ICA decomposition with and without denoising (Figure 5).
It is often suggested that ICA also has denoising properties (McKeown et al., 2003; Griffanti et al., 2014), which need be clarified. On one hand, if just a few artifactual components are removed from the signal (as in FIX cleaning), the effect of ICA is primarily to remove structured noise, and not thermal (random) noise. On the other hand, if the ICA decomposition is used to keep and examine just a few ICs that appear anatomically consistent with RSNs, then indeed most of thermal noise is also removed in that process (Beckmann and Smith, 2004). However, group-level ICA suffers from its own limitations (Cole et al., 2010) and is not necessarily the appropriate analysis tool for all studies, and seed-based analyses are expected to benefit greatly from prior denoising using MP-PCA. Another study using MP-PCA denoising prior to task fMRI analysis reported an increase of 60% in SNR and improved statistics and extent of the activation (Ades-Aron et al., 2020).
Interestingly, MP-PCA denoising alone or even combining with GSR had almost no effect in reducing intra-group variability (Figure 6, columns D, F vs. A), and it also had no contribution to the detection of between-group difference in pipelines D and F compared with the variants without it (pipelines A and B, Figure 8).
Head motion during fMRI acquisitions is one of the major confounding factors that leads to artificial correlation compromising the interpretation of rs-fMRI data (Van Dijk et al., 2012; Maknojia et al., 2019). However, compared with human studies where head motion is common, rodent studies are less impacted by this confound due to the restraint achieved by a fixation setup with ear bars and a bite bar and the use of anesthesia (Pan et al., 2015). In our datasets, no apparent head motion was observed by visual inspection of time courses except for two datasets, which were discarded. Moreover, based on the framewise displacement (FD) analysis (Power et al., 2012), head motion in our rat datasets was extremely small with displacements less than 1 μm (<0.3% of voxel size) for 75% datasets and less than 7 μm (<2% of voxel size) for all (Supplementary Figure 7), whereas in human data, it is on the order of several tenths of a millimeter (∼5–10% of voxel size) (Power et al., 2012). Motion correction was, therefore, skipped in the proposed data processing pipeline since it has been shown to introduce spurious correlations (Grootoonk et al., 2000; Chuang et al., 2019; Sirmpilatze et al., 2019).
The FIX-based artifact auto-classification has already been applied in human and mouse fMRI datasets (Salimi-Khorshidi et al., 2014; Zerbi et al., 2015). In this work, this automatic artifact removal approach was for the first time implemented for rat data. After being trained in a manually classified dataset, the classifier showed a high accuracy in identifying artifact components from rs-fMRI signal in an untouched test dataset. Nonetheless, Table 2 shows that there is a trade-off between the recall and precision in the fully automatic classification of artifact components, which means that it is not possible to achieve both very high recall and precision with one FIX threshold. However, in practice, this problem could be addressed by half automated classification in which two auto-classifications are first performed with, respectively, low and high thresholds (20 and 70 for instance), and then the difference between their two artifact lists are manually examined. In this way, by manually classifying a small portion of ICA components (∼24%), we were able to achieve a very high classification accuracy in a short time. Note that the training set was cleaned aggressively in order to give flexibility in aggressiveness/conservatism for test datasets by adjusting the threshold. This classifier is available along with the rest of the pipeline code3. The impact of ICA cleaning alone on individual datasets suggested that this approach is not sufficient for a comprehensive mitigation of artifacts in rat rs-fMRI data (Figure 7).
Although controversial, GSR is still commonly used in the analysis of rs-fMRI data (Falahpour et al., 2018) due to its capability of reducing the effects of respiration and motion on FC estimates (Birn, 2012; Yan et al., 2013; Power et al., 2014) and enhancing the spatial specificity of positive correlations (Fox et al., 2009). Here, we found that pipelines B and F that included GSR but no CL had little effect in reducing intra-group variability compared with the baseline protocol A and could also not systematically address pronounced artifacts in individual connectivity matrices. However, combined with ICA-based cleaning, GSR reduced most within-group variability in the healthy CTL group for both Cohorts 1 and 2 and revealed differences in connectivity between the CTL and STZ groups most consistently across timepoints. The benefits of CL + GSR have previously been highlighted in human fMRI data to strengthen the association between FC and behavior (Li et al., 2019).
Taken separately, minimizing within-group variability could favor pipelines that wipe out any relevant information-carrying signal while favoring widespread between-group differences could select a pipeline with the highest false-positive rate. Satisfying both criteria concomitantly though prevents either scenario. As far as low intra-group variability being potentially a signature of over-cleaning, we underline that the average FC in the control group (Figure 9, bottom left halves of each matrix) displays expected features of strong inter-hemispheric connectivity between L/R regions as well as, for instance, anti-correlation between midline regions of the DMN and the sensorimotor system, as previously reported (Gozzi and Schwarz, 2016), which argue in favor of preserving expected baseline FC with our optimized pipeline. As mentioned, MP-PCA denoising results in more clear-cut ICA decomposition where fewer components are labeled as structural noise and removed. Over-cleaning is perhaps instead manifest in the pipeline variant without MP-PCA denoising (pipeline C: ICA cleaning + GSR) that achieved the minimal intra-group variability due to more ICs being labeled as structured noise and removed (see Figure 6). The positive influence of MP-PCA denoising on downstream ICA decomposition and identification of stronger and cleaner RSNs suggests that such an approach may also be beneficial for identifying relevant variability within a control group (Bergmann et al., 2020).
A pipeline’s ability to reveal true group differences assumes that such true differences exist, which may be challenging to ascertain. First, we underline that the proposed PIRACY pipeline revealed group differences consistently across timepoints and across integration of additional datasets (Cohort 2) to the initial Cohort 1. In contrast, other pipelines either revealed no differences at all (without CL nor GSR) or differences that were inconsistent across timepoints and across dataset pooling, and thereby indeed pointed to false positives. For example, pipeline E (DN + CL) showed overwhelming group differences only at 6 weeks (Cohort 1), which disappeared in Cohort 1 + 2 analysis, and nothing otherwise. Pipeline F (DN + GSR) showed no differences at 2 and 13 weeks and overwhelming group differences at 6 and 21 weeks (Cohort 1), with only the latter being sustained in the Cohort 1 + 2 analysis.
In terms of the reliability of the differences pointed by the PIRACY pipeline, while false positives cannot be excluded due to the relatively small sample size, we stress that group differences were indeed expected at all chosen timepoints and the pattern that PIRACY showcased for these changes agreed with previous literature on this animal model assessed from histology, behavior, and volumetry (Lester-Coll et al., 2006; Kraska et al., 2012; Knezovic et al., 2015). The FC analysis resulting from Cohort 1 and its pertinence with respect to the animal model and to AD have been published separately (Tristão Pereira et al., 2021). Briefly, STZ rats in Cohort 1 exhibited not only altered FC but also intra-axonal damage and demyelination (assessed using diffusion MRI) in brain regions typical of AD, in a temporal pattern of acute injury, transient recovery/compensation, and chronic degeneration. The non-monotonic pattern in FC changes was characterized by initial hyper-connectivity and impaired network dissociation, followed by later hypo-connectivity, consistent with patterns found in pre-clinical AD and mild cognitive impairment (Dickerson et al., 2005; Schultz et al., 2017). The switch occurred between the 6- and 13-week timepoints. FC in the CTL group was consistent with previous reports that midline regions of the DMN are anti-correlated with the sensorimotor system (Gozzi and Schwarz, 2016). In STZ animals, these anti-correlations between DMN and LCN were initially reduced (2 and 6 weeks), suggesting less efficient network dissociation and brain processing (Fox et al., 2005). Concomitantly, hyperconnectivity within the DMN was found in STZ rats at these early timepoints. By 13 weeks, however, STZ rats exhibited reduced connectivity (or hypoconnectivity) in regions typically involved in AD, also consistent with eventual memory impairment in this animal model and reduced FC in AD patients (Binnewijzend et al., 2012; Franzmeier et al., 2019). While Cohort 1 had too few datasets at 21 weeks to identify group differences, the combination of Cohorts 1 + 2 in the context of the current work revealed patterns of sustained inter-hemispheric hypoconnectivity particularly between posterior brain and DMN regions at this timepoint. In the context of sustained glucose hypometabolism, these non-monotonic trends–also reported in behavioral studies of this animal model (Knezovic et al., 2015) as well as in human AD (Dickerson et al., 2005; Fortea et al., 2011, 2014; Pegueroles et al., 2017; Sierra-Marcos, 2017; Dong et al., 2020)–suggest a compensatory mechanism, possibly recruiting ketone bodies, that allows a partial and temporary repair of brain structure and function.
Limitations
First, free breathing of animals was chosen as an easier setup, best suited for longitudinal studies. While it certainly brings in higher respiratory variations to fMRI data than artificial ventilation, the breathing rate was carefully monitored and was relatively stable during the 10′ of fMRI acquisition. Breathing fluctuations were further mitigated at the level of ICA cleaning and GSR. Second, animals breathed a consistent mixture of air and oxygen, but SpO2 was not monitored. While we did maintain a strictly consistent experimental protocol throughout the study in terms of gas mixture, anesthesia levels, and timings, varying SpO2 levels might confound the BOLD signal and thereby the results. We note, however, that a recent meta-analysis has shown that, even when blood gas levels are monitored, their impact on the BOLD signal is not consistent/systematic across studies (Steiner et al., 2020), which makes it difficult to account for. Finally, how our findings generalize to completely different datasets, particularly in terms of disease model, remains to be investigated.
Conclusion
We conclude that the PIRACY processing pipeline for rat rs-fMRI data proposed herein, which includes MP-PCA denoising, a FIX auto-classification and cleaning of structured artifacts uncovered by ICA, and GSR, allowed to greatly reduce the within-group variability and improve the detection of between-group differences at the same time. This data processing pipeline, therefore, has strong potential to improve the sensitivity and reproducibility of rs-fMRI studies on rat models of disease and injury.
Data Availability Statement
The pipeline code and FIX rat classifier are available at https://github.com/Mic-map/PIRACY. The data pertaining to Cohort 1 is available at https://openneuro.org/datasets/ds003520.
Ethics Statement
The animal study was reviewed and approved by the Cantonal Service for Veterinary Affairs (Vaud, CH).
Author Contributions
IJ and RG conceived and designed the experiments. YD performed the experiments. YD and TY analyzed the data. All authors wrote the manuscript.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
The authors thank Analina Da Silva, Mario Lepore, and Stefan Mitrea for technical assistance with animal experiments and acknowledge access to the facilities and expertise of the CIBM Center for Biomedical Imaging, a Swiss Research Center of excellence founded and supported by Lausanne University Hospital (CHUV), University of Lausanne (UNIL), Ecole Polytechnique Fédérale de Lausanne (EPFL), University of Geneva (UNIGE), and Geneva University Hospitals (HUG).
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2021.602170/full#supplementary-material
Footnotes
- ^ https://fsl.fmrib.ox.ac.uk/fsl
- ^ https://www.nitrc.org/projects/whs-sd-atlas
- ^ https://github.com/Mic-map/PIRACY
References
Ades-Aron, B., Lemberskiy, G., Veraart, J., Golfinos, J., Fieremans, E., Novikov, D. S., et al. (2020). Improved task-based functional MRI language mapping in patients with brain tumors through marchenko-pastur principal component analysis denoising. Radiology 298, 365–373. doi: 10.1148/radiol.2020200822
Aquino, K. M., Fulcher, B. D., Parkes, L., Sabaroedin, K., and Fornito, A. (2019). Identifying and removing widespread signal deflections from fMRI data: rethinking the global signal regression problem. bioRxiv [Preprint]. doi: 10.1101/662726
Avants, B. B., Epstein, C. L., Grossman, M., and Gee, J. C. (2008). Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. [Special issue on the third international workshop on biomedical image registration – WBIR 2006]. Med. Image Anal. 12, 26–41. doi: 10.1016/j.media.2007.06.004
Bajic, D., Craig, M. M., Mongerson, C. R. L., Borsook, D., and Becerra, L. (2017). Identifying rodent resting-state brain networks with independent component analysis. Front. Neurosci. 11:685. doi: 10.3389/fnins.2017.00685
Beckmann, C. F., and Smith, S. M. (2004). Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23, 137–152. doi: 10.1109/TMI.2003.822821
Bergmann, E., Gofman, X., Kavushansky, A., and Kahn, I. (2020). Individual variability in functional connectivity architecture of the mouse brain. Commun. Biol. 3:738. doi: 10.1038/s42003-020-01472-5
Binnewijzend, M. A. A., Schoonheim, M. M., Sanz-Arigita, E., Wink, A. M., van der Flier, W. M., Tolboom, N., et al. (2012). Resting-state fMRI changes in Alzheimer’s disease and mild cognitive impairment. Neurobiol. Aging 33, 2018–2028. doi: 10.1016/j.neurobiolaging.2011.07.003
Birn, R. M. (2012). The role of physiological noise in resting-state functional connectivity. Neuroimage 62, 864–870. doi: 10.1016/j.neuroimage.2012.01.016
Birn, R. M., Diamond, J. B., Smith, M. A., and Bandettini, P. A. (2006). Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. Neuroimage 31, 1536–1548. doi: 10.1016/j.neuroimage.2006.02.048
Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409
Buckner, R. L., Krienen, F. M., and Yeo, B. T. T. (2013). Opportunities and limitations of intrinsic functional connectivity MRI. Nat. Neurosci. 16, 832–837. doi: 10.1038/nn.3423
Caballero-Gaudes, C., and Reynolds, R. C. (2017). Methods for cleaning the BOLD fMRI signal. Neuroimage 154, 128–149. doi: 10.1016/j.neuroimage.2016.12.018
Calhoun, V. D., Golay, X., and Pearlson, G. (2000). “Improved fMRI slice timing correction: interpolation errors and wrap around effects,” in Proceedings of the 9th Annual Meeting of ISMRM, Vol. 1, Denver, 810.
Chuang, K.-H., Lee, H.-L., Li, Z., Chang, W.-T., Nasrallah, F. A., Yeow, L. Y., et al. (2019). Evaluation of nuisance removal for functional MRI of rodent brain. Neuroimage 188, 694–709. doi: 10.1016/j.neuroimage.2018.12.048
Cole, D. M., Smith, S. M., and Beckmann, C. F. (2010). Advances and pitfalls in the analysis and interpretation of resting-state FMRI Data. Front. Syst. Neurosci. 4:8. doi: 10.3389/fnsys.2010.00008
Damoiseaux, J. S., Rombouts, S. A. R. B., Barkhof, F., Scheltens, P., Stam, C. J., Smith, S. M., et al. (2006). Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U.S.A. 103, 13848–13853. doi: 10.1073/pnas.0601417103
Dickerson, B. C., Salat, D. H., Greve, D. N., Chua, E. F., Rand-Giovannetti, E., Rentz, D. M., et al. (2005). Increased hippocampal activation in mild cognitive impairment compared to normal aging and AD. Neurology 65:404. doi: 10.1212/01.wnl.0000171450.97464.49
Does, M. D., Olesen, J. L., Harkins, K. D., Serradas-Duarte, T., Gochberg, D. F., Jespersen, S. N., et al. (2019). Evaluation of principal component analysis image denoising on multi-exponential MRI relaxometry. Magn. Reson. Med. 81, 3503–3514. doi: 10.1002/mrm.27658
Dong, J. W., Jelescu, I. O., Ades-Aron, B., Novikov, D. S., Friedman, K., Babb, J. S., et al. (2020). Diffusion MRI biomarkers of white matter microstructure vary nonmonotonically with increasing cerebral amyloid deposition. Neurobiol. Aging 89, 118–128. doi: 10.1016/j.neurobiolaging.2020.01.009
Falahpour, M., Nalci, A., and Liu, T. T. (2018). The effects of global signal regression on estimates of resting-state BOLD fMRI and EEG vigilance correlations. bioRxiv [Preprint]. doi: 10.1101/433912
Filippini, N., MacIntosh, B. J., Hough, M. G., Goodwin, G. M., Frisoni, G. B., Smith, S. M., et al. (2009). Distinct patterns of brain activity in young carriers of the APOE-ε4 allele. Proc. Natl. Acad. Sci. U.S.A. 106, 7209–7214. doi: 10.1073/pnas.0811879106
Fornito, A., and Bullmore, E. T. (2010). What can spontaneous fluctuations of the blood oxygenation-level-dependent signal tell us about psychiatric disorders? Curr. Opin. Psychiatry 23, 239–249. doi: 10.1097/YCO.0b013e328337d78d
Fortea, J., Sala-Llonch, R., Bartrés-Faz, D., Lladó, A., Solé-Padullés, C., Bosch, B., et al. (2011). Cognitively preserved subjects with transitional cerebrospinal fluid ß-amyloid 1-42 values have thicker cortex in Alzheimer’s disease vulnerable areas. Biol. Psychiatr. 70, 183–190. doi: 10.1016/j.biopsych.2011.02.017
Fortea, J., Vilaplana, E., Alcolea, D., Carmona-Iragui, M., Sánchez-Saudinos, M.-B., Sala, I., et al. (2014). Cerebrospinal fluid β-amyloid and phospho-tau biomarker interactions affecting brain structure in preclinical Alzheimer disease. Ann. Neurol. 76, 223–230. doi: 10.1002/ana.24186
Fox, M. D., and Greicius, M. (2010). Clinical applications of resting state functional connectivity. Front. Syst. Neurosci. 4:19. doi: 10.3389/fnsys.2010.00019
Fox, M. D., and Raichle, M. E. (2007). Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700–711. doi: 10.1038/nrn2201
Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. U.S.A. 102:9673. doi: 10.1073/pnas.0504136102
Fox, M. D., Zhang, D., Snyder, A. Z., and Raichle, M. E. (2009). The global signal and observed anticorrelated resting state brain networks. J. Neurophysiol. 101, 3270–3283. doi: 10.1152/jn.90777.2008
Franzmeier, N., Ren, J., Damm, A., Monté-Rubio, G., Boada, M., Ruiz, A., et al. (2019). The BDNF Val66Met SNP modulates the association between beta-amyloid and hippocampal disconnection in Alzheimer’s disease. Mol. Psychiatr. 26, 614–628. doi: 10.1038/s41380-019-0404-6
Gozzi, A., and Schwarz, A. J. (2016). Large-scale functional connectivity networks in the rodent brain. Neuroimage 127, 496–509. doi: 10.1016/j.neuroimage.2015.12.017
Grandjean, J., Schroeter, A., Batata, I., and Rudin, M. (2014). Optimization of anesthesia protocol for resting-state fMRI in mice based on differential effects of anesthetics on functional connectivity patterns. Neuroimage 102(Pt 2), 838–847. doi: 10.1016/j.neuroimage.2014.08.043
Grieb, P. (2016). Intracerebroventricular streptozotocin injections as a model of Alzheimer’s disease: in search of a relevant mechanism. Mol. Neurobiol. 53, 1741–1752. doi: 10.1007/s12035-015-9132-3
Griffanti, L., Dipasquale, O., Laganà, M. M., Nemni, R., Clerici, M., Smith, S. M., et al. (2015). Effective artifact removal in resting state fMRI data improves detection of DMN functional connectivity alteration in Alzheimer’s disease. Front. Hum. Neurosci. 9:449. doi: 10.3389/fnhum.2015.00449
Griffanti, L., Douaud, G., Bijsterbosch, J., Evangelisti, S., Alfaro-Almagro, F., Glasser, M. F., et al. (2017). Hand classification of fMRI ICA noise components. Neuroimage 154, 188–205. doi: 10.1016/j.neuroimage.2016.12.036
Griffanti, L., Salimi-Khorshidi, G., Beckmann, C. F., Auerbach, E. J., Douaud, G., Sexton, C. E., et al. (2014). ICA-based artefact removal and accelerated fMRI acquisition for improved resting state network imaging. Neuroimage 95, 232–247. doi: 10.1016/j.neuroimage.2014.03.034
Grootoonk, S., Hutton, C., Ashburner, J., Howseman, A. M., Josephs, O., Rees, G., et al. (2000). Characterization and correction of interpolation effects in the realignment of fMRI time series. Neuroimage 11, 49–57. doi: 10.1006/nimg.1999.0515
Henson, R., Büchel, C., Josephs, O., and Friston, K. (1999). The slice-timing problem in event-related fMRI. Neuroimage 9:125.
Jonckers, E., Audekerke, J. V., De Visscher, G., Van der Linden, A., and Verhoye, M. (2011). Functional connectivity fMRI of the rodent brain: comparison of functional connectivity networks in rat and mouse. PLoS One 6:e18876. doi: 10.1371/journal.pone.0018876
Kasper, L., Bollmann, S., Diaconescu, A. O., Hutton, C., Heinzle, J., Iglesias, S., et al. (2017). The physio toolbox for modeling physiological noise in fMRI data. J. Neurosci. Methods 276, 56–72. doi: 10.1016/j.jneumeth.2016.10.019
Knezovic, A., Osmanovic-Barilar, J., Curlin, M., Hof, P. R., Simic, G., Riederer, P., et al. (2015). Staging of cognitive deficits and neuropathological and ultrastructural changes in streptozotocin-induced rat model of Alzheimer’s disease. J. Neural Transm. (Vienna) 122, 577–592. doi: 10.1007/s00702-015-1394-4
Kraska, A., Santin, M. D., Dorieux, O., Joseph-Mathurin, N., Bourrin, E., Petit, F., et al. (2012). In vivo cross-sectional characterization of cerebral alterations induced by intracerebroventricular administration of streptozotocin. PLoS One 7:e46196. doi: 10.1371/journal.pone.0046196
Kruger, G., and Glover, G. H. (2001). Physiological noise in oxygenation-sensitive magnetic resonance imaging. Magn. Reson. Med. 46, 631–637. doi: 10.1002/mrm.1240
Kuang, L.-D., Lin, Q.-H., Gong, X.-F., Cong, F., Sui, J., and Calhoun, V. D. (2018). Model order effects on ICA of resting-state complex-valued fMRI data: application to schizophrenia. J. Neurosci. Methods 304, 24–38. doi: 10.1016/j.jneumeth.2018.02.013
Lester-Coll, N., Rivera, E. J., Soscia, S. J., Doiron, K., Wands, J. R., and de la Monte, S. M. (2006). Intracerebral streptozotocin model of type 3 diabetes: relevance to sporadic Alzheimer’s disease. J. Alzheimers Dis. 9, 13–33.
Li, J., Kong, R., Liégeois, R., Orban, C., Tan, Y., Sun, N., et al. (2019). Global signal regression strengthens association between resting-state functional connectivity and behavior. Neuroimage 196, 126–141. doi: 10.1016/j.neuroimage.2019.04.016
Li, Y.-O., Adalı, T., and Calhoun, V. D. (2007). Estimating the number of independent components for functional magnetic resonance imaging data. Hum. Brain Mapp. 28, 1251–1266. doi: 10.1002/hbm.20359
Liu, T. T., Nalci, A., and Falahpour, M. (2017). The global signal in fMRI: nuisance or information? Neuroimage 150, 213–229. doi: 10.1016/j.neuroimage.2017.02.036
Maknojia, S., Churchill, N. W., Schweizer, T. A., and Graham, S. J. (2019). Resting State fMRI: going through the motions. Front. Neurosci. 13:825. doi: 10.3389/fnins.2019.00825
Mayer, G., Nitsch, R., and Hoyer, S. (1990). Effects of changes in peripheral and cerebral glucose metabolism on locomotor activity, learning and memory in adult male rats. Brain Res. 532, 95–100.
McKeown, M. J., Hansen, L. K., and Sejnowski, T. J. (2003). Independent component analysis of functional MRI: what is signal and what is noise? Curr. Opin. Neurobiol. 13, 620–629.
Murphy, K., Birn, R. M., and Bandettini, P. A. (2013). Resting-state fMRI confounds and cleanup. Neuroimage 80, 349–359. doi: 10.1016/j.neuroimage.2013.04.001
Murphy, K., Birn, R. M., Handwerker, D. A., Jones, T. B., and Bandettini, P. A. (2009). The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage 44, 893–905. doi: 10.1016/j.neuroimage.2008.09.036
Murphy, K., and Fox, M. D. (2017). Towards a consensus regarding global signal regression for resting state functional connectivity MRI. Neuroimage 154, 169–173. doi: 10.1016/j.neuroimage.2016.11.052
Pan, W.-J., Billings, J. C. W., Grooms, J. K., Shakil, S., and Keilholz, S. D. (2015). Considerations for resting state functional MRI and functional connectivity studies in rodents. Front. Neurosci. 9:269. doi: 10.3389/fnins.2015.00269
Pawela, C. P., Biswal, B. B., Hudetz, A. G., Schulte, M. L., Li, R., Jones, S. R., et al. (2009). A protocol for use of medetomidine anesthesia in rats for extended studies using task-induced BOLD contrast and resting-state functional connectivity. Neuroimage 46, 1137–1147. doi: 10.1016/j.neuroimage.2009.03.004
Pegueroles, J., Vilaplana, E., Montal, V., Sampedro, F., Alcolea, D., Carmona-Iragui, M., et al. (2017). Longitudinal brain structural changes in preclinical Alzheimer’s disease. Alzheimers Dement. 13, 499–509. doi: 10.1016/j.jalz.2016.08.010
Power, J. D., Barnes, K. A., Snyder, A. Z., Schlaggar, B. L., and Petersen, S. E. (2012). Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage 59, 2142–2154. doi: 10.1016/j.neuroimage.2011.10.018
Power, J. D., Mitra, A., Laumann, T. O., Snyder, A. Z., Schlaggar, B. L., and Petersen, S. E. (2014). Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage 84, 320–341. doi: 10.1016/j.neuroimage.2013.08.048
Power, J. D., Schlaggar, B. L., and Petersen, S. E. (2015). Recent progress and outstanding issues in motion correction in resting state fMRI. Neuroimage 105, 536–551. doi: 10.1016/j.neuroimage.2014.10.044
Powers, D. (2011). Evaluation: from precision, recall and F-measure to ROC, informedness, markedness & correlation. J. Mach. Learn. Technol. 2, 37–63.
Reynaud, O., da Silva, A. R., Gruetter, R., and Jelescu, I. O. (2019). Multi-slice passband bSSFP for human and rodent fMRI at ultra-high field. J. Magn. Reson. 305, 31–40. doi: 10.1016/j.jmr.2019.05.010
Saad, Z. S., Gotts, S. J., Murphy, K., Chen, G., Jo, H. J., Martin, A., et al. (2012). Trouble at rest: how correlation patterns and group differences become distorted after global signal regression. Brain Connect. 2, 25–32. doi: 10.1089/brain.2012.0080
Salimi-Khorshidi, G., Douaud, G., Beckmann, C. F., Glasser, M. F., Griffanti, L., and Smith, S. M. (2014). Automatic denoising of functional MRI data: combining independent component analysis and hierarchical fusion of classifiers. Neuroimage 90, 449–468. doi: 10.1016/j.neuroimage.2013.11.046
Schultz, A. P., Chhatwal, J. P., Hedden, T., Mormino, E. C., Hanseeuw, B. J., Sepulcre, J., et al. (2017). Phases of hyperconnectivity and hypoconnectivity in the default mode and salience networks track with amyloid and tau in clinically normal individuals. J. Neurosci. 37, 4323–4331. doi: 10.1523/jneurosci.3263-16.2017
Shoham, S., Bejar, C., Kovalev, E., and Weinstock, M. (2003). Intracerebroventricular injection of streptozotocin causes neurotoxicity to myelin that contributes to spatial memory deficits in rats. Exp. Neurol. 184, 1043–1052. doi: 10.1016/j.expneurol.2003.08.015
Sierra-Marcos, A. (2017). Regional cerebral blood flow in mild cognitive impairment and Alzheimer’s disease measured with arterial spin labeling magnetic resonance imaging. Int. J. Alzheimers Dis. 2017:5479597. doi: 10.1155/2017/5479597
Sirmpilatze, N., Baudewig, J., and Boretius, S. (2019). Temporal stability of fMRI in medetomidine-anesthetized rats. Sci. Rep. 9:16673. doi: 10.1038/s41598-019-53144-y
Sladky, R., Friston, K. J., Tröstl, J., Cunnington, R., Moser, E., and Windischberger, C. (2011). Slice-timing effects and their correction in functional MRI. Neuroimage 58, 588–594. doi: 10.1016/j.neuroimage.2011.06.078
Smith, S. M. (2002). Fast robust automated brain extraction. Hum. Brain Mapp. 17, 143–155. doi: 10.1002/hbm.10062
Smith, S. M., Jenkinson, M., Woolrich, M. W., Beckmann, C. F., Behrens, T. E. J., Johansen-Berg, H., et al. (2004). Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23(Suppl. 1), S208–S219. doi: 10.1016/j.neuroimage.2004.07.051
Smith, S. M., Miller, K. L., Moeller, S., Xu, J., Auerbach, E. J., Woolrich, M. W., et al. (2012). Temporally-independent functional modes of spontaneous brain activity. Proc. Natl. Acad. Sci. U.S.A. 109, 3131–3136. doi: 10.1073/pnas.1121329109
Smith, S. M., Miller, K. L., Salimi-Khorshidi, G., Webster, M., Beckmann, C. F., Nichols, T. E., et al. (2011). Network modelling methods for FMRI. Neuroimage 54, 875–891. doi: 10.1016/j.neuroimage.2010.08.063
Steiner, A. R., Rousseau-Blass, F., Schroeter, A., Hartnack, S., and Bettschart-Wolfensberger, R. (2020). Systematic review: anaesthetic protocols and management as confounders in rodent blood oxygen level dependent functional magnetic resonance imaging (BOLD fMRI)–Part A: effects of changes in physiological parameters. Front. Neurosci. 14:577119. doi: 10.3389/fnins.2020.577119
Tristão Pereira, C., Diao, Y., Yin, T., da Silva, A. R., Lanz, B., Pierzchala, K., et al. (2021). Synchronous nonmonotonic changes in functional connectivity and white matter integrity in a rat model of sporadic Alzheimer’s disease. Neuroimage 225:117498. doi: 10.1016/j.neuroimage.2020.117498
Tsurugizawa, T., Djemai, B., and Zalesky, A. (2019). The impact of fasting on resting state brain networks in mice. Sci. Rep. 9:2976. doi: 10.1038/s41598-019-39851-6
Van Dijk, K. R. A., Sabuncu, M. R., and Buckner, R. L. (2012). The influence of head motion on intrinsic functional connectivity MRI. Neuroimage 59, 431–438. doi: 10.1016/j.neuroimage.2011.07.044
Veraart, J., Novikov, D. S., Christiaens, D., Ades-aron, B., Sijbers, J., and Fieremans, E. (2016). Denoising of diffusion MRI using random matrix theory. Neuroimage 142, 394–406. doi: 10.1016/j.neuroimage.2016.08.016
Wang, Y., and Li, T.-Q. (2015). Dimensionality of ICA in resting-state fMRI investigated by feature optimized classification of independent components with SVM. Front. Hum. Neurosci. 9:259. doi: 10.3389/fnhum.2015.00259
Weber, R., Ramos-Cabrer, P., Wiedermann, D., van Camp, N., and Hoehn, M. (2006). A fully noninvasive and robust experimental protocol for longitudinal fMRI studies in the rat. Neuroimage 29, 1303–1310. doi: 10.1016/j.neuroimage.2005.08.028
Yan, C.-G., Cheung, B., Kelly, C., Colcombe, S., Craddock, R. C., Di Martino, A., et al. (2013). A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics. Neuroimage 76, 183–201. doi: 10.1016/j.neuroimage.2013.03.004
Zalesky, A., Fornito, A., and Bullmore, E. T. (2010). Network-based statistic: identifying differences in brain networks. Neuroimage 53, 1197–1207. doi: 10.1016/j.neuroimage.2010.06.041
Keywords: data processing pipeline, global signal regression, ICA, resting state – fMRI, functional connectivity, rat – brain, denoising
Citation: Diao Y, Yin T, Gruetter R and Jelescu IO (2021) PIRACY: An Optimized Pipeline for Functional Connectivity Analysis in the Rat Brain. Front. Neurosci. 15:602170. doi: 10.3389/fnins.2021.602170
Received: 02 September 2020; Accepted: 26 February 2021;
Published: 26 March 2021.
Edited by:
Amir Shmuel, McGill University, CanadaReviewed by:
Wen-Ju Pan, Emory University, United StatesJoanes Grandjean, Radboud University Medical Center, Netherlands
Copyright © 2021 Diao, Yin, Gruetter and Jelescu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ileana O. Jelescu, ileana.jelescu@epfl.ch