majorly under construction and in development at the same time!

There are probably as many papers as there are preprocessing protocols, although that might be changing with the focus on reproducibility and tools that standardize pipelines like fMRIprep "A Robust Preprocessing Pipeline for fMRI Data" put out by Russell Poldrack's Lab at Stanford U and used at the Stanford Center for Reproducible Neuroscience (CRN) .

Although preprocessing is becoming a button-pushing matter, reducing the number of decision points you encounter, it is still good to understand what you're doing (how you are using your tools) and doublecheck your data -- someone even published an analysis of how the easier such complex analyses has become, the more mistakes there were in the analyses (I have to find that paper again). On the other hand, fMRI data analyses involve so many technical decisions that go beyond the expertise of the average non-methodologist systems neuroscience researcher. Acknowledging these problems most pipeline tools like fMRIprep and FSL also produce copious reports which the researcher is urged to review. 

Other tools might not have set pipelines but modules that allow you more flexibility in what is done. Here, it is imperative to keep abreast of methodological papers that raise issues that come up with more knowledge (usually simulation studies and pipeline comparisons) -- 

e.g. doing signal temporal filtering (whether band-pass or better yet high-pass -- see Shirer et al. 2015 NeuroImage) then motion regression reintroduces unwanted frequencies back into the signal (Hallquist et al. 2013 NeuroImage) -- but then so does doing it the other way around (Lindquist et al.2018 HBM)!

PREPROCESSING STEPS

Sources of noise

• participant head movements

• participant physiological noise

  • heart/blood pulsing (leads to global noise from motion and local noise from change in amount of blood and hence oxygen concentration)

  • breathing-related changes (changes in arterial CO2 concentration, movement from breathing)

• slow drifts (scanner instabilities)

Spatial realignment (within modality)

Slice timing correction/realignment

blah blah

The repetition time (TR), the time to acquire a 3D volume, is made up of 2D slices/planes. A standard TR for fMRI is around 2 sec, which means there is a difference between 2 sec in the time when the first slice was acquired and the last (assuming no parallel imaging). To compensate for the fact that slices are acquired at different points in time/make it seem like they were acquired at the same time in order to increase statistical power of models (more so for event-related than block-design experiments), we do slice timing correction/realignment. Not everyone does this...prob not needed for block design or really slow TRs but generally seems to help for event-related designs (Sladky et al., 2011).

When to do it in the pipeline? Depends on slice acquisition order and amount of movement. Prob best after spatial realignment especially if sequential slice acquisition -- so this also is a reminder that slicing timing correction might not fully fudge this so you'd have residual signal at the wrong time, e.g. when subject movement causes signal to spill across slices.

My big Q has always been: how do we know we've done it right, or rather wrong, mainly used the wrong slice timing info (or maybe wrong order in the preprocessing pipeline? Mikaël says he's seen ICA pick incorrect slice timing correction in form of striations in the ICA maps. Otherwise...??

Susceptibility distortion correction (SDC)

Huh, what's that?

MRIs are susceptible to local resonant frequencies which can cause inhomogeneities in the main static/constant homogeneous (! kind of) magnetic field (B0), which can lead to distortions (spatial mismapping) and signal loss particularly pronounced at tissue/fluid boundaries -- and this is particularly the case for echo-planar images (EPIs), which takes longer to acquire an entire 2-D slice (entire 2-D slice is acquired following a single excitation if "single-shot").

The amount of susceptibility distortions is proportional to 

• the susceptibility of an object 

  • e.g. ferromagnetic objects like metal implants or braces cause more distortions

• field strength 

  • e.g. 7T more susceptible than 3T

  • why patients with fixed metal implants are better scanned at 1.5T than 3T

• echo time (time between excitation and readout)

  • reduce SD by using shorter TE (less time for dephasing)

  • why have less SD with fast spin-echo (SE) than gradient echo (GRE -- check out Port & Pomper 2000 for minimizing SD in GRE)

• inverse to bandwidth (gradient strength and direction)

  • avoid narrow bandwidth

SDC is routinely applied for DTI/DWI analyses but not fMRI -- why oh why is that?? SDC is being recognized as a growing issue in the fMRI community but still not part of many lab protocols.

What's with the "susceptibility" in "susceptibility distortion"? Susceptibliity is actually a technical term with its own conventional Greek symbol (χ) that refers to how much something becomes magnetized when placed in a magnetic field.

SDC Approaches

fieldmap

This involves calculating the constant and linear components of the scanner's field map using info about the displacement of each point in an image -- which can be derived from the difference of phase images acquired at different echo times (TE; time between excitation and acquisition) -- in order to correct EPIs, which are based on magnitude.

Therefore, this approach usually requires the following

• images

  • two magnitude images at two different echo times ("double-echo" sequence)

  • differential phase image (difference of the phase images at the two echo times), which provides the crucial info about the displacement of each point in an image

• info

  • effective echo spacing (actual echo spacing adjusted for parallel imaging, if used, phase oversampling)

    • FSL's topup requires the readout time = echo_spacing * (matrix_lines * partial_Fourier/acceleration_factor) - 1)

      • in Siemens scanner-speak: echo spacing = "dwell time", full Fourier = partial Fourier of 1, acceleration factor = "GRAPPA"

Tools available to apply this approach:

• SPM FieldMap toolbox

• AFNI's epi_b0_correct.py

• FSL's FUGUE

• BrainVoyager's anatabacus, the document for which is peppered with wonderful details/explanations about SD and SDC à la the wonderful way of Rainer Goebel!

opposite phase-encoded images ("blip-up/blip-down")

DWIs are usually acquired with spin echo (SE; refocusing pulse to reverse the off-resonance precession, which causes a signal peak or "echo" as the spins re-align) EPI acquisition sequences, and naturally extra pairs of SE images in opposite phase-encoding directions are collected to estimate field inhomogeneities leading to such distortions in order to correct them. fMRI is usually acquired using gradient echo (GRE; pulses used to vary resonance frequence spatially in an object, usually linearly in the x/y/z planes) EPI acquisition sequences -- does that mean that GRE instead of SE images should be used for SDC? Check out my ISMRM (International Society for Magnetic Resonance in Medicine) paper on whether to use GRE or SE blip pairs to SDC fMRIs.

Tools available

• SPM ACID toolbox: HySCO routine

• FSL topup

• Computational Morphometry Toolkit (CMTK) (also see the slides below for a quick-and-dirty intro on how to use it)

"fieldmap-less" estimation

Estimate field distortion map using diffeomorphic (nonlinear) registration SyN. This approach used to be included in fmriprep, but branched off as a separate Python library, SDCFlows.

Coregistration (between image modalities)

Different modalities have different intensities for different regions (T1 - dark CSF, light WM but T2 - light CSF, dark WM) so minimizing squared differences, usually based on intensities, as with realignment within a modality is inappropriate -- use mutual information (2D histogram of mutual information between two images -- the sharper the histogram is along the diagonal, the better the alignment)

Removing physiological noise

RETROICOR noise modeling

TAPAS Physio Toolbox (Kasper et al. 2017 JoN Meth)

CompCor (Behzadi et al. 2007) - PCA of noise in WM and CSF

tICA: PESTICA (Beall & Lowe 2007)

sICA: CORSICA (Perlbarg et al. 2007)

sICA: ICA-AROMA (Pruim et al. 2015)

FIX-ICA (Salimi-Khorshidi et al. 2014) included in FSL, requires classifier training

multi-echo acquisitions + ME-ICA (Kundu et al. 2012)

Caballero-Gaudes et al. (2017 NeuroImage) - Methods for cleaning the BOLD fMRI signal

external recordings (pulsometer, breathing belts)

Removing motion-related artifacts and other noise

Note that many of the references below focus on connectivity analyses of resting-state fMRI data, but they should also apply to task-based fMRI data analyses.

• global signal regression (GSR): whole brain

  • assumption: global signal dominated mostly by noise, whether due to physiological noise, participant movement, or the scanner

  • apparently quite controversial -- improves most noise-reducing strategies (Parkes et al. 2018) but can ...

    • introduce artifactual anti-correlations because reduces centers correlations between voxels around zero (Fox et al. 2009 J Neurophysiol)

    • include signal tightly coupled GM voxels

    • cause artifactual group differences (Saad et al. 2012 Brain Connect)

  • alternative to GSR: Aquino et al. (2020, NI)

• regressing out timeseries in WM and CSF

• realignment (head motion) parameters (RP or HMP)

• signal spikes -- check out the 2018 neurostars discussion about these spike strategies among some of the major minds behind the major MRI analyses tools including Bob Cox of AFNI and Cyril Pernet of SPM.

  • strategies

    • motion scrubbing (also sometimes referred to as motion censoring) : remove volumes where frame displacement (FD; i.e. volume-by-volume movement -- different ways of calculating this: see Yan et al. 2013 NI) exceeds some threshold (Power et al. 2012, 2014)

      • advantages:

      • disadvantages: 

    • despiking: truncate amplitude of spikes

      • advantages:

        • retains data

        • preserves time points/does not mess with time-course structure

        • voxel-by-voxel

      • disadvantages: bad data still there, just limited

    • spike regression: model each spike as a separate regressor (and hence separate delta function) in your regression model (Lemieux et al. 2007 MRI; Satterthwaite et al. 2013 NI)

      • equivalent to deleting the time point from the data -- but this does not mean the results are the same (see Ciric et al. 2017)

      • advantages:

      • disadvantages: cannot be done voxelwise

    • CompCor: temporal PCA estimation of noise in time courses in WM and CSF (Behzadi et al. 2017 NI extract components of temporal variance of a single mask covering WM+CSF; Muschelli et al. 2014 NI run PCA separately on WM and CSF) 

      • advantages: no data censoring

      • disadvantages:

    • ICA-AROMA: ICA (FSL's melodic) on smoothed data with a priori criteria to detect noise components (Pruim et al. 2015)

      • advantages:

      • disadvantages:

      • Ciric et al. (2017) and Burgess et al. (2016) both find that ICA-based approaches need to be complemented by GSR to capture motion artifacts.

    • .... the list continues

  • modeling with external physiological recordings

    • advantages: 

    • disadvantages:

      • does not effectively capture global noise

  • impact of these strategies affect (see Parkes et al. 2018 for evaluation of different strategies and measures)

    • test-retest reliability

    • QC-FC distance-dependence: movement can artifactually increase close-range functional connection strengths and reduce long-distance connections (Power et al. 2012; Satterthwaite et al. 2012)

Ciric et al. (2017 NI) and Parkes et al. (2018 NI) compare common strategies, offer guidelines on how to decide what strategy, and how to assess impact on data quality. Parkes et al. give a procedure. Ciric et al. find that for their dataset and measures, overall the best overall for their dataset and the measures they looked at (with lots of caveats) was (1) ICA-AROMA + GSR or alternatively (2) GSR + RP + WM + CSF + spike censoring

ICA-AROMA: ICA-based Automatic Removal Of Motion Artifacts

ICA-AROMA (Pruim et al. 2015 NeuroImage, available as standalone Python code on GitHub and included in fMRIprep) identifies and removes motion-related ICA components from FMRI data. It uses a small set of theoretically motivated features, specifically meant for head-related motion, which avoids having to do any training of the noise classifier that is applied to your dataset. 

You can include the the results of ICA-AROMA in what they call the "non-aggressive" method, which involves using the EPI BOLD images that have been directly noise-corrected (in fMRIprep these have the affix desc-smoothAROMAnonaggr_bold -- note they are already spatially smoothed which is required by AROMA unlike the non-AROMA files (smoothing is not part of fmriprep) and note also that ICA-AROMA only outputs 2mm isotropic images because it only uses the MIN152 2mm template, as of time of writing 2021) , or you can use your non-AROMA-changed EPI BOLD images and do noise regression modeling (including the principle components identified as regressors in your GLM) which they refer to as the "aggressive" method. Either way, it is imperative to review what components have been identifed as noise to make sure that real signal has not been filtered out. 

One diagnostic that might indicate that ICA-AROMA has taken out real signal is to compare the tSNR ( = mean(signal)/SD(signal) ) in GM tissue of your original EPIs with the AROMA'd EPIs -- if tSNR is much lower for the AROMA'd EPIs then check your AROMA components more carefully.

For each participant, fMRIprep creates a summary HTML report in which there is a section, "ICA Components classified by AROMA,"  when ICA-AROMA is used, listing each ICA component localized on a glass brain (maximum intensity projection -- look for bright red or deep blue blobs) with the time series of the component (top line plot) and the frequency of the of the signal identified (bottom line plot), and classified as noise (red) or "real" signal (green). 

Griffanti et al. (2017 NeuroImage) put out a guide to identify ICA noise components in fMRI. This thread on the NeuroStars forum which started off about comparing different noise removal models with ICA-AROMA results and ensuing problems (e.g. report that ICA-AROMA is not very adapted to multiband images)

But here are some things to look out for (THANK YOU, MIKAËL NAVEAU for walking me through this! But any mistakes below are mine). Check that:

• the components identified as noise (red line plots) are indeed noise and do not account for an inordinate amount of the signal

  • normal "noise"

    • slow linear trends or signal drift due to scanner instabilities (Smith et al. 1999 NeuroImage): 0.0 - 0.015 Hz and spatially distributed

    • physiological noise: frequency line plots should show (a) peak(s), indicating regularity, which is what you would associate with breathing or heart beats

      • breathing 

        • normal rate ranges -- 6 years old and older : slow 4.5 - 6.5 bpm (0.075 - 0.11 Hz), normal 12 - 20 bpm (0.20 - 0.33 Hz)

        • can also lead to low frequency fluctuations (<0.1 Hz) because change in arterial CO2

      • cardiac pulsation

        • average adult normal heart rate ranges: 50 bpm (0.83 Hz), 100 bpm (1.67 Hz)

        • localized near large arteries/veins (e.g. sagittal sinus, middle cerebral artery), edges of the brain, lateral ventricals, sulci

        • By the way, that adds up to about 2.5 billion times that the heart beats in an average lifetime (The Physics Factbook)!

    • tissue pulsation or CSF flow due to heart pulses -- localized around huge arteries and ventricles

    • movement: e.g. clusters around the edge of the brain or stripey looking clusters (usually interaction of head motion with scanner artifacts)

• the components identified as real data (green line plots) are indeed plausible real data

  • actual BOLD signal usually has a regular and low frequency but can go up to 0.01-0.1 Hz and clusters should be large and not too many

  • depending on your task, e.g. for visual tasks involving a motor response, look for components capturing activity in the visual cortex and in the motor cortex

• check that what got tagged as noise (red) has been effectively removed in the components tagged as "good" (green) 

  • e.g. if a couple of components has captured slow scanner drifts, ideally this drift should not appear in any of the signal components

• no severe striation effects

  • e.g. if the slice order information for slice timing correction was incorrect, ICA can pick this up and show this as lines spaced at regular intervals, corresponding to the slices that were specified in the wrong order -- however, this can sometimes also be a multiband artifact

If any components are incorrectly classified or do not capture the full noise, you can tweak the ICA-AROMA parameters and rerun and iteratively come to a refined set of components for each individual -- or you can model the components like Mikaël does -- either way, a time-consuming process.

Examples

Real-data examples! Although ICA-AROMA was optimized for head-related movements, and it will more likely classify a component as noise if it correlates with known motion components, it will also pick up other noises like low-frequency scanner drifts and physiological noise. Here are some examples (please send me a message if I am wrong or you have any comments/feedback!)

Coregistration

Normalization

Spatially warp all the apples into oranges so we can compare them all like oranges.

Smoothing

Handling the "noise" introduced by subjects moving (but not suddent movements/spikes in signal), imprecise normalization and coregistration, etc. by assuming that neighboring voxels have similar effects, by spatially blurring the image using a moving Gaussian to "stabilize" the signal in each voxel. Basically adjusting the intensity values in each voxel factoring in the distance-weighted intensities of its neighbors (the closer the higher the weight/taken into consideration more), and the size of the Gaussian is referred to by the full-width half-maximum (FWHM), the width/extent at half the maximum height of the peak.

What FWHM to choose? The bigger the FWHM, the blurrier your image, i.e. spatially inaccurate but the stronger your statistical power. Rule of thumb: 2x or 3x the size of your image resolution (so for 2x2x2mm images, go for around 4-6-mm FWHM). But this should also depend on the size of the regions you are interested in...you don't necessarily want a big FWHM that blurs across different anatomical regions. 

Online resources

• Wikibooks on neuroimaging data processing