This talk will describe the main steps of functional MRI data preprocessing based on the blood-oxygenation level dependent (BOLD) contrast. The numerous methods available for each step, and corresponding parameter selection, causes that the amount of possible preprocessing workflows can be enormous, which may lead to substantial variability in the quality of the preprocessed data and final results. We will establish some simple guidelines for adequate preprocessing since there is no ‘optimal’ preprocessing pipeline, but there are incorrectly applied methods, emphasizing that the workflow must be decided according to the characteristics of each dataset and the research question. We will briefly introduce several platforms that can help researchers to design the preprocessing pipeline, automatize its execution, and facilitate data quality assessment. These tools can foster reproducibility, and ensure transparent reporting of methodological details.
The amount of possible data preprocessing workflows in fMRI data analyses can be enormous, which may lead to substantial variability in the quality of the preprocessed data and conclusions from fMRI results (Carp, 2012; Churchill et al., 2015). Since there is no ‘optimal’ preprocessing pipeline, but there are incorrectly applied methods, researchers must be aware of the pros and cons of each decision and understand that the elements of the preprocessing pipeline and their relative order must be decided according to the characteristics of each dataset (e.g. conventional single-echo fMRI, multi-echo fMRI) and the research question (e.g. whole-brain parcellation-based analysis vs. mesoscopic fMRI)
The preprocessing pipeline could start by voxelwise despiking (e.g. 3dDespike in AFNI or the ArtRepair toolbox) to improve the accuracy of volume registration (Jo et al., 2013), even though it can also result in less apparent head motion in the images (Power et al., 2017). Next, a block of several operations, including slice-timing correction, volume registration, correction of magnetic field distortions and physiological noise correction, could be applied. The relative order of these operations is controversial since they are all interconnected to each other. The majority of fMRI data is collected using 2D multi-slice acquisition sequences where one slice is acquired at a time, either in an interleaved or a sequential order. Therefore, different parts of the 3D image are acquired at different times, up to TR*(1-1/Nslices) seconds. Slice timing correction aims to compensate this offset by interpolating the voxel time series as if all of the slices were acquired at the time of a reference slice. The effectiveness of slice timing correction is still debated. Some researchers indicate that its use is not detrimental (Sladky et al., 2011), whereas others suggest it is not necessary with fast acquisitions (TR ≤ 1 sec) or not effective with long TR (TR ≥ 2.5 s).
Realignment or volume registration compensates the misalignment between images of the time series due to head motion by realigning all images to a reference image, either a single image or the mean of the time series. This step is also known as motion correction, although it is important to understand that the effects of motion are not actually corrected in this step. Realignment involves an initial estimation of motion, usually assuming a rigid body transformation (i.e. the position of head can only change in 3 translational (x, y and z) and rotational (pitch, roll and yaw) directions) that are obtained by optimizing a cost function (e.g. weighted least-squares or normalized correlation ratio are usually employed), and optionally creating the realigned images through spatial interpolation. If processing time is not stringent and spatial interpolation is to be applied at this stage of the preprocessing, higher-order interpolation methods (e.g. heptic, sinc, spline, Fourier) should be preferred over more simple interpolations (linear, cubic).
Slice timing and realignment interact closely and depending on their order. If slice-timing correction is performed first, the amount of estimated motion can be smaller than the actual one due to the temporal interpolation (Power et al., 2017) motion-related changes can propagate across time. Furthermore, any actual through-plane motion will result in a mismatch between the nominal and actual timing of data acquisition. On the other hand, if volume realignment is performed first including spatial interpolation, data acquired at one time point will be combined with data from a different slice, and thus the nominal acquisition time may not again match the actual acquisition timing, reducing the effectiveness of slice timing correction. For these reasons, the order of slice acquisition, either sequential or interleaved, and the existence of slice gap also matter. For instance, when acquiring slices in an interleaved fashion, the time between excitation of adjacent slices is one half TR. Therefore, there is more opportunity for motion to affect the slice timing and introduce artefacts. A possible solution is to combine both volume registration and slice-timing correction in one step (Bannister et al., 2007; Roche et al., 2011), although this approach is rarely applied.
The most common fMRI sequence, gradient-echo echo-planar imaging (EPI), suffers from artefacts in regions close to tissue-air interfaces due to magnetic field inhomogeneities. These artefacts manifest as a reduction in the intensity of the fMRI signal and geometric distortions along the phase encoding direction that is used to sample the k-space. The signal dropout can be reduced by acquiring data at shorter echo times, with smaller voxels, or with multi-echo fMRI sequences. Geometric distortions can be corrected based on magnetic field maps (Jezzard and Balaban, 1995; Hutton et al., 2002) or two acquisitions with reversed phase-encoding directions (Andersson et al., 2003). These images must be acquired close in time to the actual functional run to correct. Notably, head motion also interacts with the correction of susceptibility-induced distortions because the effects of head motion in regions with susceptibility artefacts do not obey a rigid-body model and voxels will be warped differently and dynamically. Consequently, combining the spatial transformations of correction for geometric distortions and volume registration in a single step is also recommendable. Even so, the interaction of volume registration with other steps becomes more relevant as multichannel receiver coils have become standard in MRI systems.
Incorporating physiological denoising into the preprocessing workflow further complicates the choice of the order of these initial steps. Whereas some authors have found that the best pipeline is to first perform volume registration, then physiological denoising and finally slice timing (Jones et al., 2008), other authors have recommended applying physiological denoising (e.g. RETROICOR (Glover et al., 2000)) prior to any step so as to accurately account for phasic changes in the fMRI signal associated with cardiac pulsations and respiratory cycles (Jo et al., 2013). This controversy depends on type of acquisition and the patterns of motion artefacts in the data. For example, in an axial interleaved acquisition with large through-plane motions due to head nodding, spatial and temporal interpolations associated with realignment and slice-timing correction will likely mix voxel time series acquired one half-TR apart. In that case, the phase of the cardiac and respiratory cycles at the time of acquisition of these voxels might be substantially different. An alternative approach is integrating the effects of motion correction into physiological denoising (Jones et al., 2008).
Coregistration
between the subject’s anatomical image to the functional data is usually performed
with 12 parameters known as a linear or affine transformation (i.e.
translation, rotation, scaling and shearing) computed based on a
between-modality optimization (e.g. mutual information), which can be thus easily
inverted. The spatial transformation is usually computed between from the reference
volume in realignment (after correcting for geometric distortions if applicable)
to the anatomical image, and then applied it to the entire functional data. Transforming the anatomical image rather than
the functional volumes is more recommended for analyses in the subject’s space.
The resulting anatomical-to-functional transformation can also be applied to extra
images resulting from tissue-based segmentations or brain parcellations to
compute tissue-based regressors for denoising (see below) or constrain subsequent
analyses within voxels of interest, for instance in grey matter. As noted
above, functional data acquired with EPI is prone to geometric distortions that
complicate the alignment with the anatomical image if the distortions are not
previously corrected. Therefore, a good approach is to acquire an
inversion-recovery (IR)-EPI volume with T1-weighted contrast with the same or
higher spatial resolution as the functional data. This T1-weighted EPI image will
have the geometric distortions matched to the functional data and thus can
serve as an intermediate reference image to improve anatomical-functional
registration (Renvall et al., 2016).
Often, the data acquired from multiple subjects is combined to generalize findings across a given population, or compare the results of a single individual with those of a population. Since individual brains are highly variable in size and shape, a common step of the preprocessing is the spatial transformation of the functional data into a standard or stereotactic coordinate system (Tailarach or MNI). Spatial normalization can be done on the functional data prior to statistical analysis (recommended in SPM and AFNI), or in the result images after statistical analyses (used in FSL). Spatial normalization is usually performed by a multistep volume-based registration to a template image (e.g. MNI305 or MNI152). First, the subject’s high-resolution anatomical image is moved to match an image template in stereotactic coordinates, which can be done through an initial affine transformation followed by a subsequent nonlinear transformation to improve the match between the anatomical image and the template. Nonlinear transformations are often described in terms of basis functions, such as the discrete cosine transform or splines used in SPM, or even as the concatenation of polynomial transformations applied in local patches of the image used in AFNI. Next, these transformations are applied to the functional data or the statistical result images. A good practice is to concatenate all of the spatial transformations (geometric distortion, volume realignment, coregistration to anatomical image, and spatial normalization to template) to create a single spatial transformation from the fMRI native space to the stereotactic space, rather than transforming and interpolating the images with each step. High-order interpolation approaches should be preferred for the interpolation of the functional data to the standard coordinate system. Visualization of the transformed images is recommended to avoid warp that are clearly not anatomically reasonable. Note that spatial normalization is not a mandatory step in the preprocessing, and can be avoided when the analyses demand high spatial accuracy or consider the variability in the shape of the brain across individuals.
The final steps of the preprocessing comprise spatial smoothing, and optionally temporal filtering or denoising. The order of spatial smoothing and denoising can be exchanged as long as the nuisance regressors used for denoising are defined from anatomical masks before spatial smoothing to avoid mixing contributions from different tissue types (Jo et al., 2013). Typically, spatial smoothing is performed with an isotropic volumetric kernel (e.g. a Gaussian kernel with a full-width half maximum larger than the voxel size), or can be limited to a given mask (e.g. 3dBlurInMask in AFNI) to avoid leaking from other tissues or non-brain voxels. In a more elaborated manner, surfaced-based spatial normalization and smoothing approaches are becoming increasingly popular as an alternative to volume-based approaches (Robinson et al., 2014, Glasser et al., 2016, Jo et al., 2007; 2008, Blazejewska et al, 2019) to diminish blurring across regions with different functionality, offering a better area localization compared with traditional volume-based methods (Coalson et al., 2018).
On the other hand, the number of existing algorithms for denoising fMRI data can be overwhelming, where each method is more or less tailored for a given type of noise. Generally speaking, most of the methods are based on the idea of nuisance regression, where a set of regressors aiming to describe artefactual or confounding components of the signal is fit to the voxel time series, and then their explained variance is removed from the data (Caballero-Gaudes and Reynolds, 2017). It is important to apply nuisance regression all at once (Lindquist et al., 2019; Jo et al., 2013; Hallquist et al, 2013), and also account for the degrees of freedom in the analysis.