Core preprocessing pipelineThe preprocessing of functional magnetic resonance imaging (fMRI) scans prior to analysis consists of, at minimum, the anatomical alignment of scans to a common space, head realignment to correct for motion, and the correction of susceptibility distortions arising from the echo-planar imaging (EPI) acquisition of functional scans. The core preprocessing pipeline in Rodent Automated Bold Improvement of EPI Sequences (RABIES) carries each of these steps with state-of-the-art processing tools and techniques (detail of each step in Supplementary Table 1).To conduct common space alignment, structural images, which should be acquired along the EPI scans, are initially corrected for inhomogeneities and then registered together to allow the alignment of different MRI acquisitions. This registration is conducted by generating an unbiased data-driven template through the iterative linear and non-linear registration of each image to the dataset consensus average (using ANTs’ SyN registration algorithm40 with Mattes Mutual Information and ANTs cross-correlation objective functions for linear and non-linear stages, respectively), where the average gets updated at each iteration to provide an increasingly representative dataset template (https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction)41. The finalised template after the last iteration provides a representative alignment of each MRI session to a template that shares the acquisition properties of the dataset (e.g., brain shape, field of view, anatomical contrast,…), making it a stable registration target for cross-subject alignment. This newly-generated unbiased template is then itself registered to an external reference atlas to provide both an anatomical segmentation and a common space comparable across studies defined from the provided reference atlas.The remaining preprocessing involves the EPI images. A volumetric EPI image is first derived using a trimmed mean across the EPI frames, after an initial motion realignment step. Using this volumetric EPI as a target, the head motion parameters are estimated by realigning each EPI frame to the target using a rigid registration. To correct for EPI susceptibility distortions, the volumetric EPI is first subjected to an inhomogeneity correction step and then registered non-linearly to the anatomical scan from the same MRI session, which allows the calculation of the required geometrical transforms for recovering brain anatomy42. This data-driven approach to distortion correction was opted for as a default technique to maximise accessibility, like fMRIPrep, since it does not require additional acquisitions which may not be available (e.g., B0 field map or reversed phase encode image). Finally, after calculating the transformations required to correct for head motion and susceptibility distortions, both transforms are concatenated into a single resampling operation (avoiding multiple resampling) which is applied at each EPI frame, generating the preprocessed EPI timeseries in native space7. Preprocessed timeseries in common space are also generated by further concatenating the transforms allowing resampling to the reference atlas.The preprocessing workflow of RABIES is illustrated in Fig. 1a, and each associated module is described in detail in the RABIES documentation (https://rabies.readthedocs.io/en/stable/). When structural images are not provided by the dataset, the volumetric EPI replaces the structural image during the generation of the unbiased template.Datasets for preprocessing benchmarkingThe robustness of RABIES preprocessing was evaluated by sampling a range of datasets from different sites to attempt to capture representative data samples across the rodent fMRI community. We accessed 20 mouse fMRI datasets (totalising 318 C57Bl/6 J mice including both males and females), as well as 3 independent rat fMRI datasets (totalising 166 rats including Wistar (both male and females), Long-Evans (male only), and Sprague-Dawley (male and females) strains). No sex-specific analyses were conducted since the information in this regard was incomplete in the original publication of the datasets. All datasets are openly accessible, and acquisition details for each dataset can be found in their respective original publication (see Supplementary Table 2). All data was acquired with approval from local ethical authorities.For mouse preprocessing using structural scans, we used the DSURQE ex-vivo MRI atlas43,44,45,46, and the Fischer 344 in vivo atlas47 was used for rat preprocessing. 17 mouse datasets did not include structural scans (see Supplementary Table 2), allowing for testing the workflow without structural scans, and for which we generated a mouse EPI atlas for the preprocessing of mouse datasets lacking structural images (see Methods section 3 below).Generation of a EPI atlasWe selected 5 mouse datasets (i.e., 117_Cryo_mediso_v, 7_Cryo_med_f1, 7_Cryo_med_f2, 7_Cryo_aw_f, 7_RT_med_f in Supplementary Table 2) from the multi-site study by Grandjean and colleagues9. We selected those datasets for their complete brain coverage, minimal susceptibility distortions and variable anatomical contrast. Using the RABIES workflow, each EPI scan was first corrected for inhomogeneity and an unbiased template was generated combining all 5 datasets. The unbiased template was then registered to the DSURQE mouse atlas using a non-linear registration40 to propagate brain masks and anatomical labels onto the newly generated EPI template. This EPI atlas is openly available online (https://doi.org/10.5281/zenodo.5118030), and used by default by RABIES as reference atlas when executing the EPI-only preprocessing workflow.Preprocessing quality controlDuring preprocessing, a set of images is generated to conduct visual quality control of each failure-prone preprocessing step (Supplementary Fig. 1). To visualise the intensity inhomogeneity correction of either the structural or the EPI image, the input MRI image is shown before and after the correction together with the overlap of the brain mask used for the final correction, thus allowing to visualise the quality of inhomogeneity correction. Failures can be noticed either if important inhomogeneities remain or if brain masking fails. Then for the remaining registration steps, the report shows the overlap between the moving image and the target image. These registration steps include alignment to the generated unbiased template, the alignment between the unbiased template and the reference common space atlas, and the registration between anatomical and functional images for susceptibility distortion correction.The entire quality control report from each dataset was visually inspected to record the number of failures at each registration step. Registration failures were considered when there was a clear misalignment of at least one brain structure. Minimal misalignments that subsisted because of EPI distortions and poor anatomical contrast were tolerated, given that few EPIs have sufficient contrast to achieve specific anatomical alignment and complete correction for distortions through non-linear registration. In Supplementary Table 2, additional notes on preprocessing quality are provided to detail dataset-specific limitations (e.g., partial misalignment to the reference atlas in highly distorted regions). All reports inspected are provided in this paper in the supplementary files.Confound correction workflowThe workflow for confound correction is structured following best practices found in human literature, largely based on recommendations from Power and colleagues27. Here, we detail the execution of the confound correction workflow, which requires specific considerations regarding the order and combination of multiple correction steps. These considerations aim to avoid inefficient corrections, or the re-introduction of confounds into the data at later steps28. Several options are provided, and, if selected, are applied in the following sequence:
1.
Frame censoring: Frame censoring temporal masks are derived from framewise displacement and/or DVARS thresholds, and applied first on both BOLD timeseries before any other correction step to exclude signal spikes which may bias downstream corrections, in particular, detrending, frequency filtering, and confound regression27.
a.
Censoring with framewise displacement: Apply frame censoring based on a framewise displacement threshold. The frames that exceed the given threshold, together with 1 back and 2 forward frames, will be masked out17.
b.
Censoring with DVARS: The DVARS time course is z-scored (i.e., subtract the temporal mean and divide by the temporal standard deviation), and frames with values above a threshold of 2.5 (i.e., above 2.5 standard deviations from the mean) are removed. z-scoring and outlier detection is repeated within the remaining frames, iteratively, until no more outliers are detected, to obtain a final set of frames post-censoring. This second censoring strategy was implemented to address residual spiking artefacts which were not well detected with a one-shot framewise displacement censoring approach (Supplementary Fig. 17).
2.
Detrending: Linear (or quadratic) trends are removed from the timeseries. Detrended timeseries are obtained by performing ordinary-least square (OLS) linear regression, where regressors consist of time together with an intercept (time-squared is also included as a regressor if removing quadratic trends).
3.
ICA-AROMA: Motion sources can be cleaned using a rodent-adapted ICA-AROMA classifier48 (Methods section 6). ICA-AROMA is applied prior to frequency filtering to remove further motion effects that can result in ringing after filtering26,48.
4.
Frequency filtering: For highpass and/or lowpass filtering of timeseries, censored timepoints are first simulated while presenting the frequency profile of the timeseries to allow for the application of a Butterworth filter.
a.
Simulating censored time points: frequency filtering requires particular considerations when applied after frame censoring since conventional filters cannot handle missing data (censoring results in missing time points). To address this issue, we implemented a method described in Power et al.27, allowing the simulation of data points while preserving the frequency composition of the data. This method relies on an adaptation of the Lomb-Scargle periodogram, which allows estimating the frequency composition of the timeseries despite missing data points, and from that estimation, missing time points can be simulated while preserving the frequency profile49.
b.
Butterworth filter: Following the simulation, frequency filtering is applied using a 3rd-order Butterworth filter, and timepoints can be removed at each edge of the timeseries to account for edge artefacts following filtering27. Edge artefacts were visualised in simulated data (Supplementary Fig. 18), and from those observations, we conclude that removing 30 s at each edge and using a 3rd order filter is optimal for a highpass filter at 0.01 Hz. After frequency filtering, the temporal mask from censoring is re-applied to remove simulated time points.
5.
Confound regression: For each voxel timeseries, a selected set of nuisance regressors (Supplementary Table 5) are modelled using OLS linear regression and their modelled contribution to the signal is removed.$$\begin{array}{c}{{\mathbf{\beta }}}={{\rm{OLS}}}({{\bf{X}}},{{\bf{Y}}})\\ {{{\bf{Y}}}}_{{{\bf{CR}}}}={{\bf{X}}}{{\mathbf{\beta }}}\\ \hat{{{\bf{Y}}}}={{\bf{Y}}}-\,{{{\bf{Y}}}}_{{{\bf{CR}}}}\end{array}$$
(1)
Regressed timeseries \(\hat{{{\bf{Y}}}}\) are obtained with Eq. (1) where \({{\bf{Y}}}\) is the timeseries, \({{\bf{X}}}\) is the set of nuisance timecourses and \({{{\bf{Y}}}}_{{{\bf{CR}}}}\) is the confound timeseries predicted from the model at each voxel (\({{{\bf{Y}}}}_{{{\bf{CR}}}}\) is a time by voxel 2D matrix). Prior to the regression, the nuisance timecourses were subjected to the same frame censoring, detrending and frequency filtering, which were applied to the BOLD timeseries to avoid the re-introduction of previously corrected confounds27,28.
6.
Intensity scaling: Signal amplitude must be scaled to allow for comparison of network amplitude estimates between scans32,50. The following options are available within RABIES:
a.
Grand mean: Timeseries are divided by the mean intensity across the brain, and then multiplied by 100 to obtain percent BOLD deviations from the mean. The mean intensity of each voxel is derived from the linear coefficient (β) from the intercept computed during detrending.
b.
Voxelwise mean: Same as grand mean, but each voxel is independently scaled by its own mean signal.
c.
Global standard deviation: Timeseries are divided by the total standard deviation across all voxel timeseries.
d.
Voxelwise z-scoring: Each voxel is divided by its standard deviation (i.e., z-scoring).
7.
Smoothing: Timeseries are spatially smoothed using a Gaussian smoothing filter.
Importantly, each confound correction step (with the exception of linear detrending) is optional when using RABIES to allow for adapting the correction strategy to specific dataset needs.Rodent-adapted ICA-AROMAThe original code for the algorithm (https://github.com/maartenmennes/ICA-AROMA) was adapted to function without the hard-coded human priors for anatomical masking and parameter thresholds for component classification. Following an initial independent component analysis (ICA) decomposition of the data using FSL’s MELODIC algorithm51, four features are extracted from each ICA component spatial map for classification: (1) the temporal correlation between the component and motion parameters, (2) the fraction of the frequency spectrum within the high-frequency range, and the fraction of the component within (3) the cerebrospinal fluid (CSF) mask and (4) the brain edge mask. The component is classified as motion if the CSF or high-frequency content fractions are above a given threshold, or if classified by a pre-trained linear classifier combining the brain edge fraction and motion correlation. To adapt the original algorithm with RABIES, the CSF mask is inherited from the rodent reference atlas and the edge mask is automatically generated from the brain mask, the threshold for high frequency content was increased as rodent can express higher BOLD frequencies, particularly under medetomidine52, and the linear classifier was retrained. To select the new parameters, we manually classified motion and network components derived from a set of scans from the REST-AWK group anaesthetised under a medetomidine-isoflurane mixture, and selected parameters to successfully classify clear motion components while avoiding false classification of brain networks.Standardisation of variance voxelwiseObserving the distribution of signal variability across the brain (i.e., BOLD variability) allows one to identify confound signatures if present, whereas relatively uncorrupted scans present a predominantly homogeneous distribution of signal variability. To attempt mitigating the impact of confounds observed from the BOLD variability map (Fig. 2), we implemented a correction strategy to enforce a homogeneous variance distribution across voxels. To do so, timeseries are first scaled voxelwise by their standard deviation (yielding standardised variance), and then timeseries are re-scaled to preserve the original total standard deviation of the entire 4D timeseries (i.e., the global standard deviation does not change). This is to avoid changing the overall scale of image intensity, which is handled separately during confound correction (see Methods section 5). Importantly, this is different from z-scoring, where timeseries are standardised to unit variance. Enforcing standardised variance distribution effectively downscales the contribution of voxels with higher variance, which typically corresponds to a confound signature.Connectivity analysis implemented within RABIESHere we detail how each connectivity analysis available within RABIES can be conducted following the completion of the confound correction workflow.
Seed-based connectivity29: An anatomical seed of interest is selected, and the mean time course within the seed is correlated (Pearson’s r) with each voxel timeseries to measure the brain connectivity map for that seed. In RABIES, this is repeated for each scan independently, to obtain a subject-level connectivity map that can then be used for hypothesis testing.
Whole-brain connectivity matrix30: the whole brain parcellation during RABIES preprocessing is used to extract a mean time course at each parcel. Then, all parcel timecourses are cross-correlated (Pearson’s r) to obtain a whole-brain matrix describing the connectivity strength between every region pair. As with SBC, in RABIES, a connectivity matrix is computed for each scan.
Group independent component analysis (ICA)51: the timeseries from all scans are concatenated in common space, and ICA decomposition is conducted on those concatenated timeseries to identify spatial sources of covariance (components) prevalent in the entire group. Among those components, resting-state networks as well as sources of confounds can be identified (example in Supplementary Fig. 19). RABIES uses FSL’s MELODIC ICA algorithm51.
Dual regression31,32: Dual regression is a technique building on the group ICA algorithm which allows the modelling of previously detected group-level components back onto single scans. The dual regression algorithm consists of two consecutive linear regression steps. First, components from group ICA are regressed against an individual scan timeseries to obtain an associated time course for each component. Then, a second regression using these component timecourses is regressed again against the scan’s timeseries, allowing to finally derive scan-specific spatial maps for each corresponding ICA component. Those linear coefficients can allow for describing individual-specific features of resting-state networks, such as amplitude or anatomical shape. To include network amplitude, the timeseries of the first regression are variance-normalised prior to the second regression, to weight the variance explained by each component into the spatial maps32. This method is applied to all dual regression analyses in this paper.
Datasets included for data quality characterisation and guidelinesTo explore signatures of confounds in mouse fMRI, we first relied on a dataset with a group of anaesthetised mice under a medetomidine-isoflurane mixture equipped with mechanical ventilation (i.e., controlling for motion) (N = 10 females) and a group of physically-restrained, freely-breathing, awake mice (N = 9 females)53. The group of awake mice was conditioned to the restraint system for 1 week prior to scanning. fMRI (gradient echo EPI, 0.18 × 0.15 mm and 90 × 60 matrix in-plane, 28 0.4 mm slices, TR = 1.2 s, 180 volumes) were acquired on an 11.7 T Bruker animal scanner. For each subject, data was acquired over 2 sessions where 3 consecutive fMRI scans were acquired per session, giving a total of 60 and 53 functional scans for the anaesthetised and awake group, respectively (one awake scan was lost during registration quality control). Furthermore, to benchmark the quality control metrics and guidelines across a representative range of data quality found with mouse fMRI, we further leveraged a multicenter dataset previously harmonised by Grandjean and colleagues (17 acquisition sites, N = 15 scans per site, including both males and females)9. All datasets were preprocessed using only EPI scans, as the multicenter study did not provide structural scans.Group-ICA priors for dual regressionThe group-ICA components used as priors for dual regression (available online https://zenodo.org/record/5118030/files/melodic_IC.nii.gz) were obtained using FSL’s MELODIC algorithm51. MELODIC was run with 30 components on the combined data from the anaesthetised and awake mouse groups of the REST-AWK dataset, where confound correction consisted of highpass filtering at 0.01 Hz, FD and DVARS censoring, and confound regression was applied with the 6 rigid motion parameters together with white matter (WM)/CSF mask signals. The group-ICA maps were visually inspected, and inspired by suggestions from Zerbi and colleagues19, each component was classified into a rodent brain network, confound source or other (Supplementary Fig. 19). More specifically, confound components display signatures of motion (i.e., edge effects), bilateral asymmetry, high loadings into WM and CSF tissues, wide-spread effects across the brain or affect only single slices. We opted for a conservative categorisation of components into confounds by selecting those which primarily displayed anatomically recognisable confound features without network-like features. For network analyses in this study, we selected two components corresponding to canonical mouse networks: (1) the somatomotor network, primarily expressed along the somatosensory and motor cortices9,19, and (2) the default mode network, primarily expressed along the cingulate and retrosplenial cortices9,19.Confound correction strategies investigatedTo assess how confound correction can improve the quality of network analysis, the REST-AWK dataset and the 17 datasets from the multi-site study9 were processed over a range of confound correction strategies. The datasets were initially corrected by applying a framewise displacement threshold of 0.05 mm and regressing the 6 motion parameters to first establish a baseline with minimal correction.We then evaluated the impact of the following additional corrections: voxelwise variance standardisation (Methods section 7), DVARS censoring, 24 motion parameters regression25, WM and CSF signal regression, aCompCor regression of the first 5 principal components, global signal regression, highpass filtering and ICA-AROMA. Lowpass filtering was not considered in this study, as lowpass inflates temporal autocorrelations54 and impacts the measures of confound correlations at the scan level. Lowpass can be used in conjunction with the proposed framework in future studies but would require re-evaluating an appropriate threshold separately.In all conditions, we applied linear detrending, grand mean scaling and spatial smoothing using a 0.3 mm Full-Width at Half Maximum Gaussian filter. When highpass filtering was applied, 30 s were removed at each end of the timeseries to avoid edge artefacts. After censoring time points, scans which had less than \(2/3\) of the original number of frames were removed from further analysis.Connectivity analyses and definition of canonical network mapsIn this study, we consider the analysis of two different mouse canonical networks, namely the somatomotor and default mode networks9,19,55, and each network was measured with either dual regression or seed-based connectivity. For seed-based connectivity analysis, the somatomotor and default mode networks were mapped using a seed in the right primary somatosensory cortex and anterior cingulate cortex, respectively. For dual regression, both networks were found among the group ICA components used as priors (Methods section 10).To conduct the proposed analysis quality control framework, we needed to define a reference brain map for each canonical network and each analysis technique. For dual regression, a representative map of each network was already defined from the group ICA priors. However, for seed-based connectivity analysis (Supplementary Figs. 15, 16), we needed to compute a representative average brain map for each network from the set of scans across datasets. To do so, we selected a subset of the datasets with various acquisition parameters (for the somatomotor network: 7_Cryo_aw_f, 7_Cryo_mediso_v, 7_RT_halo_v, 94_RT_iso_v, REST-AWK mediso; for the default mode network: 7_Cryo_mediso_v, 7_RT_halo_v, 94_RT_iso_v, 117_Cryo_mediso_v, REST-AWK mediso) and performed well on quality metrics with dual regression analysis, when using the initial confound correction strategy (framewise displacement censoring + 6 motion parameters regression). The consensus network map was obtained by taking the median connectivity across scans. The selection of datasets was repeated independently for each network to select datasets succeeding quality control for the relevant network. The resulting consensus network maps were highly similar to the corresponding ICA components.Dataset distribution reportTo complement the group statistical report during analysis quality control (Fig. 4), as well as the selection of scan-level thresholds for network specificity and confound correlation (Fig. 3), an additional figure is generated to display the sample distribution across quality metrics (Supplementary Fig. 7). This report is automatically generated along the spatiotemporal diagnosis and group statistical reports and allows visualising the relationship of both network specificity and amplitude with each measure of confound, i.e., the within scan temporal correlation with confounds, mean framewise displacement, the variance explained from confound regression and temporal degrees of freedom. For dual regression only, the overall network amplitude is summarised for each scan by taking the L2-norm across the spatial map of network amplitude.This report can be used to assess the proportion of scans which pass quality control thresholds, as well as evaluate the impact of modifying confound correction strategy on network specificity and temporal correlation with confounds. In addition, visualising the relationship between network amplitude and the three measures of confounds from the group statistical report allows for inspecting the correlational structure between connectivity and confounds. In particular, it can be useful for detecting outliers, which are automatically identified using a threshold of 3.5 on the modified Z-score56 for each measure independently (network amplitude/shape, and each confound measure). Z-scores are computed on the subset of scans which passed scan-level quality control thresholds and after removing outliers in network amplitude (Methods section 14).Outlier removal based on network amplitudeScans presenting outlier values in overall network amplitude from the distribution plot were automatically removed, as this likely indicates spurious connectivity (e.g., Supplementary Fig. 8). After removing scans which did not pass quality control, outliers were automatically detected and removed using the modified z-score, with a threshold of 3.556. If outliers were detected and removed, z-scores were re-computed on the remaining scans, and new outliers were removed. This process was repeated iteratively until no more outliers were detected to obtain a group distribution without outliers.Confound correction optimisation protocolWe designed a stepwise protocol for optimising confound correction on a per-dataset basis, which refers to a table relating data quality features to the corresponding correction strategy (Supplementary Fig. 14a). For all analyses, we selected the following scan-level quality control thresholds: minimum of 0.4 Dice overlap for network specificity and maximum of 0.25 correlation for confound effects (Fig. 3). The procedure applied for each dataset was as follow:
1.
The dataset is initially corrected only using framewise displacement censoring and regression of 6 motion parameters, and the initial data quality reports are generated.
2.
Evaluation of data quality reports:
a.
Using the spatiotemporal diagnosis, the 4 main data quality markers described in Fig. 2 are inspected across scans to identify potential issues.
b.
The dataset distribution report is inspected to assess how many scans pass the quality control thresholds and identify outliers. If fewer than 8 scans passed, the group statistical report was not considered.
c.
Inspect the group statistical report to notice absent or spurious signatures in the network variability map and evaluate the extent of confound correlation (the average correlation within the network is saved in an attached CSV file).
3.
An appropriate CR strategy is selected based on the observations from the data quality reports, following the order of priority and instructions in Supplementary Fig. 14a.
4.
One additional correction is applied at a time to evaluate the impact of individual corrections. The quality outcomes are re-evaluated as in (2), and the correction is kept only if improving quality outcomes.
5.
Repeat (3), (4) until no quality issues are left, or CR options are exhausted.
When evaluating the quality reports, the somatomotor – dual regression analysis was prioritised, but other analyses were also considered for improvements. The optimised set of corrections defined for each dataset is documented in Supplementary Fig. 14b.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.