fMRIPrep for a robust preprocessing pipeline

When dealing with neuroimaging studies, researchers often create their own preprocessing workflows that suit the individual labs and datasets. The existence of a variety of MRI analytical tools does not help with the standardization of the analysis. Recently, a group at Stanford University led by Russ Poldrack et al. came up with a robust and analysis-agnostic tool for robust and reproducible preprocessing called fMRIPrep. The software tool integrates existing software packages such as FSL, ANTs, AFNI, Nypipe, and SPM. It also works on the BIDS dataset, another useful framework to standardize neuroimaging data collection following the spirit of Open Science.

Preprocessing anatomical images

1. Correction of intensity non-uniformity with N4BiasFieldCorrection (ANTs).
2. Then skull-stripped the image with antsBrainExtraction.sh script (ANTs).
3. Run recon-all (FreeSurfer) to reconstruct 2D cortical surface from the 3D T1-weighted image.
4. Cortical grey matter segmentation is done using FreeSurfer.
5. Do non-linear spatial normalization to the standard template* using antsRegistration (ANTs).
6. Brain tissues: CSF, white matter/WM, and grey matter/GM are extracted using FAST (FSL).

* standard space follows the ICMB 152 nonlinear asymmetrical template (MNI, v2009c). You can also choose the slightly older version, MNI152 asymmetrical template 6th generation. 

Preprocessing functional images

1. For each BOLD run found in the folder, we first find its reference volume, which is usually the 3D volume at midpoint (half-way through the run); and it should be skull-stripped also.
2. Head-motion parameters are estimated using MCFLIRT (FSL).
3. Slice timing correction is performed using 3dTshift (AFNI) if slice-time information is available.
4. Fieldmap distortion correction is then applied for SDC around the air space.
5. Co-registration into the subject's T1 image is done using FLIRT with BBS and 6 d.o.f (FSL).
6. If surface reconstruction is selected, then bbresgiter (FreeSurfer) is applied for co-registration. The authors reported that this method yields the best results owing to the accuracy of GM/WM surfaces driving the process.
7. A series of transforms of BOLD image up to the co-registration can be customized using antsApplyTransforms and Lanczos interpolation.
8. BOLD images can be normalized to the MNI template of a different resolution: 1mm or 2mm, and to a different template version.
   
9. Lastly, if ICA-AROMA* is enabled, then the pipeline will perform ICA analysis and do some denoising. At this stage, preprocessed data up to Step-7 has to be "smoothened" first before feeding it into FSL-MELODIC. 

* By default, fmriPrep will produce a non-aggresive denoised and preprocessed BOLD image in the standard space at the end of ICA-AROMA. The authors suggest that non-aggresive method of denoising is recommended as it removes less signal components that share variance with the nuisance regressors.  

Extraction of nuisance regressors

By default, there is no temporal denoising done in fMRIPrep but it produces a vast variety of nuisance regressors. This allows researchers to be more flexible in their denoising strategies.

1. Component-based correction (CompCor) provides physiological noise regressors.
2. PCA on high-pass filtered data is performed to yield tCompCor (top 5% of the variable voxels in the subcortical masks without any GM components) and aCompCor (intersection between mask and CSF/WM).
3. Three nuisance regressors from CSF, WM, and whole brain signals.
4. Framewise displacements and their first derivatives to signify head movements are also obtained (DVARS).
5. Components produced by ICA-AROMA are also included.

Additional Notes on fMRI

TR and TR in MRI

When deciding a suitable MRI sequence, you have to decide the choice between short or long repetition time (TR) and echo time (TE). Different combinations of these two values produce different types of MR images. See below. To recap, the T2* image has transverse magnetization which decays much faster than would be predicted by natural atomic and molecular mechanisms. In practice, T2* can be considered as an effective/observed T2 image, whereas the T2 image is considered the natural or true T2 (i.e. T2* is always less than or equal to T2). The T2* images result principally from inhomogeneities in the main magnetic field.

Statistical power in MRI

Statistical power is defined as the probability of rejecting the null hypothesis when it is false. In most cases, it is dependent on three main factors: (1) the desired effect size; (2) alpha value or Type I error; (3) the sample size tested. Of these 3 values, the sample size is the one that is within the control of the experimenter (Desmond & Glover, 2002). Unlike in usual behavioural tasks, fMRI data of each voxel is typically normalized as percent signal % change between the Task and Rest conditions. Variability in the fMRI data consists of the within-scan and between-scan variability. Increasing scan time and the number of subjects will potentially reduce these two kinds of variation. 

The concept of statistical power in fMRI data is complex because it involves different expected activation patterns, analysis pipelines, and also scanning parameters. Loosely, some rule of thumb to obtain a decent signal, alpha = 0.05 and 80% power. For example:

a)  Total N = 12 - 20 subjects per condition.

b)  Scan time per subject = 10 - 15 min.

c)  Scan volume (time points) = 400 - 500.

Non-parametric Approach to MRI

Like the usual inferential statistics, we want to make a conclusion in our neuroimaging data that apply to the population. As the name implies, non-parametric tests do not depend on the nature of the distribution of the fMRI data (e.g. normality assumption of the residual component). Through resampling, we construct an empirical distribution from the data itself that can accurately estimate the probability of a particular observance. This is the fundamental idea of non-parametric statistics. Some of these approaches have been clearly discussed in Nichols & Holmes, 2001. The use of non-parametric statistics if the data is not normal has been discussed briefly in this blog previously. A more thorough review article that includes the history and logic behind non-parametric tests is by Winkler A, et al. (2014) in the NeuroImage journal.

Basic Concepts
To see why this makes sense, the task-based block design paradigm will be used, i.e. a subject performs a task in one block and rest in the other. If a specific location of the brain is used during tasks, we should expect a task-related activation during the task, but not during the rest period. Here, we reshuffle or scramble the data across TRs. If a particular voxel is involved in the task, we should not expect any strong correlation between the predicted BOLD response and the actual activation because the 'reshuffled data' is just random noise. Reshuffling does not change the mean and variance of the dataset, that is, we don't tamper with the data.

In analyzing task-based fMRI data, for example, you form a design matrix according to your behavioural task and conduct a GLM (using FEAT, let say) to obtain parametric estimates or PEs of the corresponding regressors. To correct for the familywise error, the multiple comparison problem across voxels, we then use GRF theory by defining cluster-forming threshold (Z = 2.30 in FEAT). According to Echlund et al, GRF performance in controlling Type-I Error is poor. In fact, it is recommended to use non-parametric tests such as permutation testing to control for Type-I Error.

Using randomise on my dataset
For my iMac machine purchased in 2011 with Intel i5 Quad-core, 3.1 GHz, it takes me 8-9 minutes to analyze a decent task-based, 100 x 100 x 63 x 250, preprocessed 4D NIFTI file using 1000x permutation. With FEAT-glm using the same GLM design matrix, it takes me 2-3 minutes.

Registration with ANTs

When I first embarked in the fMRI project, I didn't learn about other software package. Recently, I managed to compare the performance of T1 -> MNI normalization done with 3 different popular software package: ANTs, FSL, and AFNI. The performance was measured using "similarity metric" that I took from Nipype, a python-based neuroimaging package. My data came from MRI scans of stroke patients, as part of my project at the Jewish General Hospital. As a reference image, I took the MNI-152 with 1mm resolution. The left panel below shows what happened when I used MNI-152 2mm template instead. In both FSL_1 and FSL_2, I used the 2 mm template to obtain the non-linear transformation matrix. However, in FSL_1, I explicitly specified the 1mm template when I applied the non-linear warping (with applywarp) process.

There are two main types of cost function: intra-modal (least squares and normalised correlation) and inter-modal (correlation ratio and mutual information-based options).


The left panel shows the registration accuracy of the subject's T1 image to the standard MNI 1mm template of 3 different software packages. The right panel also shows what happened when the MNI 2mm template was used instead.

As a side note, the MNI template has been widely used in various neuroimaging software packages. It was historically produced in a report for the International Consortium for Brain Mapping (ICBM) back in 2001. The atlas was based on linear mapping of 152 healthy adults ranging from 18-44 years old, mostly white people. It has gone through 2 major revisions lately to increase the accuracy: the first one was in 2006 that made use of a non-linear mapping, they call this the MNI-152 6th generation. FSL uses this template. Another revision was the most recent one in 2009.

Recent development in functional MRI

More recent contributions in functional MRI have been based on graph theory (Bullmore & Sporns, 2009), which underlines the importance of seeing the brain as complex networks. This analysis of MRI data has been widely used in both healthy and dysfunctional brain networks. 

It is important for understanding normal brain development and functions, the networks involved in affect regulation, motor control and execution, and learning and memory (e.g. motor learning study, see Sami and Miall, 2013). In addition, the graph theory has been beneficial to understand the implication resulting from e.g. Alzheimer's disease, schizophrenia, and stroke. For example, Buckner et al. found a correlation between the site of the targeted regions and the location of major hubs in Alzheimer’s disease (2009). He et al. (2007) demonstrated how the frontoparietal network is implicated in the spatial neglect of stroke patients.

In order to establish a brain network using modern graph theory, there are a number of steps to be taken: 

- define the network nodes (usually from resting-state fMRI), 

- estimate association/correlation between nodes, 

- compile pairwise associations between nodes and generate an association matrix, 

- and, calculate the network characteristics. 

Yet, despite this emerging literature, many researchers using graph theoretical approaches for fMRI data have not closely examined whether their data violate the assumptions of graph-theory. One important assumption that has to be met is the assumption of stationarity. Time series or time-varying data is said to fulfill the stationarity principle if the statistical properties of a time series are time-invariant, meaning, the mean and variance do not change across time. Generally speaking, stationarity implies that the statistic parameter of interest does not change over time. Such assumption is also essential for analyzing fMRI time series in the frequency domain, as the Fourier transform is suitable for stationarity. For more than a decade, most scientists have always assumed that fMRI signals are stationary, however, one may refer to some recent works that warrant doing fMRI analyses with caution (e.g. Muhei-aldin et al., 2014; Ou et al, 2014).

Notes on Resting-state fMRI Analyses (Part II)

Fieldmap Correction and Coregistration
What is fieldmap correction?
Apart from having a poor resolution, functional images acquired using common EPI sequences suffer distortion due to magnetic field (B₀) inhomogeneity introduced by different tissue types in our heads. Such a thing occurs due to the existence of non-homogeneity in RF receive and transmit of the head coils. The more channels you have, the more inhomogeneity the image may have. The most severe inhomogeneity includes the air-bone or air-brain tissue interfaces in the sinuses in the inferior frontal gyrus and medial temporal lobes. This poses a serious effect on our data, a geometrical distortion and signal loss as depicted in Fig-1.

Magnetic field inhomogeneity can be measured with fieldmap images; which can give us a geometric distortion and signal loss. These values can then be used to compensate for the loss by geometrically unwarping the EPI images, and applying cost-function masking in registrations to ignore areas of signal loss. The correction is most useful during image co-registration as it dramatically improves the registration accuracy. Areas where signal loss has occurred unfortunately cannot be restored with any form of post-processing. In other words, it is impossible to recover time-series data in those locations.

There is no separate sequence for acquiring the fieldmap and different scanners give different images. The sequence can be EPI, Spin-echo, or Gradient-echo sequences, but it isn't recommended to use the EPI-based sequence since it will suffer the same problem. There exist 2 different methods of acquiring fieldmap images for the purpose of correction. 

When you do the fieldmap acquisition, you usually acquire two different images: a pair of magnitude images captured with different echo times, and a phase difference image (Fig-1). The acquisition can also be controlled either in the AP (j+) or PA (j-) direction. These images should be acquired in the same orientation as the target EPIs. The phase difference between the two images is proportional to the difference in echo time (ΔTE) and the B0 inhomogeneity observed. The fieldmap is calculated by taking the difference between the two-phase images, and dividing that by the echo time difference.

Method-2 is called the blip-up blip-down method, which calculates the fieldmap based on the difference in distortion between the two consecutive acquisitions. This method acquires two diffusion-weighted images (DWI) with opposite phase encoding directions, that is, the AP and AP directions. It is assumed that there is no change in the magnetic field and sudden motion during the two acquisitions. You can use TOPUP in FSL to help you do fieldmap processing using this method.

Fig-1: Images obtained from the scanner (left) are converted to get a fieldmap image (right). Red circles show distorted regions that require correction. This is a standard procedure of double gradient-echo performed in Siemens 3T scanner.

How to process this in FSL?
At the MNI, our brain imaging center uses Siemens 3T scanner, which is a good thing as FSL provides a ready-to-use tool, fsl_prepare_fieldmap, to obtain a fieldmap phase image in rad/sec. The magnitude image resembles a lower resolution version of the T1 structural (anatomical) image. FSL FUGUE, which is incorporated in FEAT, helps us to do distortion correction using this method. Both the complete and skull-stripped versions of the magnitude image and the processed phase image (rad/sec) should be defined in FEAT. FSL will then attempt to unwarp the distorted EPI image before mapping it to the structural image. The unwarp direction has to be specified and is typically given by the scanner operator depending on how the fieldmap acquisition is set.

In FSL, fieldmap correction is incorporated as part of the registration (preprocessing) pipeline. The highly accurate functional-to-structural coregistration is also called boundary-based registration or BBR  (Greve and Fischl, 2009). The method is based on changes in the intensity along the white matter boundaries instead of the less reliable grey matter boundaries. This means that an accurate segmentation of the structural image is required and bias-field correction reliably improves the accuracy. Performing BBR registration without a fieldmap correction doesn't give many benefits than the usual 6DOF method with FLIRT (Fig 2-3).

Also, there must be some grey-white intensity contrast in the EPI, though it doesn't have to be good enough for segmentation. The FSL website said since only intensities near the white-matter boundary are used by BBR, it is likely to be more robust to a range of pathologies and artefacts in the EPI or the structural.

Fig-2: Comparison of 3 situations using BBR coregistration with fieldmap correction: when there full magnitude image with the skull wasn't supplied to FSL(left); when the correct magnitude image was used but the unwarping direction was the opposite (middle); the correct BBR-registered image with a superior accuracy (right).

Fig-3: Performing coregistration of a functional image to a structural image using BBR is superior than the usual linear 6DOF registration in FSL. Note that the asterisk ( * ) sign indicates the region with severe signal loss. Without using the fieldmap correction, the corpus callosum mapping becomes inaccurate as denoted by a hex sign (#).

 

Slice Timing Correction?

Scientists more or less agree that the slice timing correction is important.  For the more recent multiband sequence, some experts said that slice timing misalignment may not have a huge impact on the analysis. In the earlier version of the Siemens WIP, I was told that the slice timing information contained in the DICOM files was not correct. This can be retrieved easily with a Matlab function. As a result, I didn't perform this correction in my fMRI paper (MB3, TR=1690 msec).

Until recently, one can deduce the slice timing information based on the CMRR Multiband protocol here. For comparison, I have included the effect of slice timing correction to my resting state data with an MB 3x acceleration measured on a single voxel.

Fig-4: Time series with and without the slice timing correction measured on a single voxel @ MNI coordinate (67,41,49).

Adding additional EVs to GLM

Additional regressors (EVs) can be added to the GLM in FSL FEAT. I find that the GUI is a bit tricky, better write a script for that. First, we have to recreate the design matrix by adding the extra regressors, using either:

       ⁍ Pointing to a file for each regressor by constructing a full model design    

       ⁍ Creating a space-delimited text-file comprising all confound EVs

Fig-5: If you click the "Full model setup", a new GUI will appear as shown on the left. Select an appropriate setting (number 1-3). You don't have to perform another temporal filtering. The temporal derivative is optional too. Another way is to construct a text file and select "Add additional confound EVs" (number 4). 

Refer to the Fig-5 above. If you click "Full model setup", a new GUI will appear. Choose the input file as 1-entry per volume (see number 1), no need to convolve it with the HRF anymore (number 2). Also, you do not have to perform another temporal filtering if the regressors are derived from the prefiltered data (number 3). This method is time-consuming. A better option is to use the additional confound EVs (number 4), just that the text file has to be space-delimited, not a comma-separated file! Using a wrong delimiter will cause FEAT to ignore these additional regressors! 

It is also safer to write a code or script rather than getting restricted with the GUI features. To create a design matrix, use the feat_model command. To manually perform GLM according to the design matrix with prewhitening, use the film_gls command. The command is equipped with sophisticated estimations of autocorrelation and Tukey tapering, which is very important for making statistical inferences in task-based fMRI.

One last note about the design matrix is about Orthogonalization. Most EVs generally are almost orthogonal, so enforcing orthogonality is not gonna help much in the results.

Denoising nuisance components with ICA

As mentioned often, rs-fMRI has one major drawback: the data is recorded at rest so it is prone to noise or artefacts. This is so because we don't have any reference pattern as we do when we perform task-based functional imaging. Hence. proper cleanup is paramount to getting a correct deduction or conclusion. This is even more serious for me who is doing learning-related rs-fMRI. One of the things I'm struggling with is choosing the best cleaning method! Technically for my project, I plan to try out the ICA denoising method since the work of our previous postdoc used a different technique.

As mentioned in the previous post, the subject-level data cleaning steps cover the following:

1. First, you perform ICA, e.g. using FSL MELODIC and identify the nuisance components. From my experience, the tool performs a pretty good job. Let the algorithm choose the best number of ICs.
2. Identify the nuisance components and run fsl_regfilt script in FSL, producing a so-called clean or denoised fMRI dataset. 
3. Now, to assess its performance after cleanup, you can either conduct another round of ICA on the residual image or compute the temporal standard deviation of the GM regions (or a specific ROI in the motor cortex).

The script above effectively removes nuisance components identified from the original dataset. If I conducted another ICA on the denoised data, the resulting components are much cleaner!

How does the script differ from the usual GLM function in FEAT? They are not quite the same. The core function of the GLM in FEAT is film_gls. The fsl_regfilt, on the other hand, uses a simple time-varying GLM function called fsl_glm, which does not perform any sophisticated modelling of temporal autocorrelation and whitening. I finally found this difference after struggling for so long!    

Fig-6: Demeaned time-series of the same voxel obtained from FEAT and ICA denoising tool.

The output of regression is called the residual image, res4d.nii.gz, which is supposed to be as clean as the denoised image, but having the mean removed (just to add back). See Fig-6 for representative time-series of a voxel from the residual outputs of each FEAT and fsl_regfilt. There are some minor differences. I think pre-whitening is not necessary because the nuisance EVs are all spatially independent from the MELODIC. 

On the other hand, I finally found out that the residual output of the fsl_regfilt is in fact the same as the one produced by the fsl_glm. The data has to be demeaned, and the design matrix regressors des_norm has to be normalized into unit variance. The FSL gurus in their forum claimed that both methods use the same GLM methods basically, but I'm not sure which GLM it was. Knowing this similarity is essential because now I can compare different methods of denoising, e.g. using WM/CSF average time-series.

To regress out nuisance components following either ICA or other noise modelling tool (e.g. RETROICOR), you can just use the simple a GLM function. Sophisticated estimation of temporal autocorrelation becomes important when you want to do statistical inference of neural activity or connectivity.

Notes on Resting-state fMRI Analyses (Part I)

{ I have done some practice on the neuroimaging analyses and I took down some interesting findings as "notes" }

Bandpass Filter, BPF
Bandpass filtering is essential to most standard resting-state (rs-fMRI) pipelines. When Biswal (1995) wrote the first rs-fMRI study, low-pass filtering with a 0.08 Hz cutoff was applied to the dataset. Biswal took a seed or ROI in the cortical motor area and conduct a whole-brain temporal correlation. He also observed that the power spectrum of the spontaneous fluctuation resides mainly < 0.1 Hz. Other earlier works, e.g. Lowe et al., 1998 and Cordes et al., 2000 follow the same footsteps.

A study focused on this question showed that only frequencies below 0.1 Hz contribute to regionally specific BOLD correlations, with faster frequencies relating to cardiac or respiratory factors. Based on this finding, the majority of spontaneous BOLD studies low-pass filter data at a cut-off of 0.08 or 0.1 Hz. [Fox & Raichle, 2007]

What is a bandpass filter?
A bandpass filter only allows certain frequency range to pass through while attenuating the frequency components outside that range. There are two cutoff frequencies associated with a bandpass filter, low and high cutoff. Roughly speaking, it can be achieved by a combination of low-pass and high-pass filtering. Some studies prefer to use the word 'low-pass' filtering on the resting-state data. This is valid because high-pass filter has already been applied to remove drift during preprocessing steps. The use of Butterworth filter is common. In Matlab, it is best to use forward and reverse algorithms to prevent phase shift and distortion. Other filter includes Gaussian filter (FSL, e.g. in Damoiseaux, et al., 2006).

Why and when should we filter?
First reason is the definition of the rs-fMRI as the low-frequency spontaneous signal fluctuation. Second reason is to avoid the influence of physiological noise, the primary non-neuronal signals that corrupt our dataset. Seed-based and ICA-based methods are the most common methods in rs-fMRI analysis. In the ICA-based method, the higher frequency artefacts will be identified as individual components. Unlike ICA, the seed-based or ROI method requires us to remove unwanted artifacts manually. There is one warning. Bandpass filtering works best only at very fast acquisition or very low TR. For a TR = 0.50 msec, for example, the highest frequency range is up to 2 Hz. In other words, aliasing is prevented up to ~ 1 Hz. This is high enough to capture both respiratory and cardiac components. For higher TR, a combination of filtering and multiple regression is used to clean up the data. For more review, see: DP Auer (2008)

What is the cutoff frequency?
The lower cutoff frequency is normally 0.01 Hz, a usual nominal value to remove drift. The earliest and most common high cutoff value is either 0.08 or 0.1 Hz. Another popular value of 0.15 Hz is used by other studies (e.g. Greicius, et al., 2003 and Ellen. et al., 2008) and it is recommended by some (e.g. Urs Braun, et al., 2012). It seems that the high cutoff value is more ambiguous, but in order to be comfortable in choosing that value, we should study the frequency content of the resting brain. I like one of the earliest studies about this topic by Cordes et al., 2001. Using a fast TR (high speed acquisition to reduce aliasing), frequency components of 0 - 0.1 Hz are found to dominate three different areas in the brain: motor, auditory, and visual. More recently, a higher value of 0.25 Hz is popular with a much faster acquisition (e.g. Boubela, et al., 2013 and Kalcher, et al., 2014).
How do we extract frequency spectrum in rs-fMRI dataset? DYI = One can use Matlab to get the time series of a voxel then use FFT to get the frequency spectrum. Alternatively, one can use ICA to separate the components and identify the components of interest, say, associated with the primary motor cortex.

Identifying Nuisance Components by ICA
Spatial ICA, which separates rs-fMRI data into spatially independent patterns of activity, has been popular as an exploratory fMRI data analysis. What does ICA give you? The spatial and temporal components, plus the frequency content of that component (McKeown et al., 2003). This can be extremely useful to check whether a certain IC is resting-state component or just noise.

Single-subject ICA Approach
There are two ways you can use ICA with your resting-state data. First, by Group ICA (gICA) to perform model-free multivariate exploratory analysis on the multiple subject dataset. This can be done through temporal concatenation. Second, to a single-subject dataset. Usually, the noise components identified here will be used in the subsequent pipeline in the GLM, that is, to regress them out from the dataset. This can be achieved in FSL by using the fsl_regfilt command after performing MELODIC.

What is the optimum # IC?
Often, the number of calculated # IC is also known as model order or ICA dimensionality. It can be freely selected up to n − 1, where n is the number of time points or volumes. It is common to set 20 - 30 using a standard pulse sequence with long TR, ~150 volumes (e.g. Calhoun, et al., 2001; Greicius, et al., 2007) or 60 (e.g., van de Ven, et al. 2007). Using a probabilistic approach, FSL-MELODIC sets the automatic dimensionality option by default so I leave this setting as it is. With this option, different participant dataset yields to a different number of components, (e.g. de Luca, et al., 2005). Under-estimating the model order may prevent us from capturing the full spectra of noise in the data. This results in a less effective method to clean the data. Also, the ICs may contain both signals and noise, making the judgment difficult. Overestimating the # IC, according to Li YO. et al., 2007, reduces the stability of the IC estimates and degrades the estimation of task-related brain activations. Another nice discussion on this topic is in a study by Abou-Elseoud et al., 2010.

How can we identify noise/nuisance ICs?
Manual visual inspection is the most direct approach. A nice publication with figures by Kelly RE Jr., et al., 2010 tells us generic rules to identify nuisance components from the ICA results. Some important points to label the components as noise:
    More than 75% of the frequency bands are > 0.1 Hz. More useful by using faster TR.
    Activation around the perimeter is highly likely due to movements.
    Activation around the ventricles, especially the lateral and fourth ventricles are usually obvious.
    Sporadic little activation blobs in the white matter.
    Activation around the major arteries or sinuses.
    Sudden spike (movement) or signal drops (artifacts).

    The separation is often not that clear-cut. If doubtful, do not remove that component.
Once these components have been identified, we can place them in GLM as regressors. Efforts are made to create an automated software to classify signals and noise and to denoise resting-state data, e.g. SOCK and FSL-FIX.

How many # ICs should I remove?
Again there is no consensus but it assumes that the more noise you throw, the cleaner your data will be. Unfortunately, to some datasets, ICA does not do the job well. This is particularly true with high TR, causing aliasing in the temporal dataset. In other words, some components contain both signal and noise, or spatially it looks noisy but temporally its frequency spectrum is around 0 - 0.15 Hz. Based on my personal experience with regular and multi-band accelerated pulse sequences (MB = 3), I often found that almost 50% of the components are noise. The challenge would be to carefully select the noise while at the same time avoiding a false positive.

Group Analysis in Functional MRI - a brief overview

Previously, I mentioned subject-level or first-level GLM analysis that is performed on individual subject data. This step gives us within-subject contrasts across parameter estimates, e.g. estimates of brain activity due to visual stimulation vs. rest. In this blog entry, a second-level analysis will be discussed which provides group-level inference on whether they are significantly different from zero, e.g., group status, behavioral performance, etc. Statistically speaking, the goal of the standard group analyses is to do inference of a set of acquired sample to a larger population. How can we perform group analysis in neuroimaging data? There are two ways we can model the "effects" of a behavioral study on the brain, that is, fixed effects and random effects. A statistic that takes both variations into account is called mixed-effects modeling in psychology or behavioral sciences. 

Consider the following scenario. You want to know the distribution of hair length between adult men and women, but doing this for all world population is just impossible. You start to select four men and women. For every single person, there is bound to be differences in hair length which represents a within-subject variability. In addition, different people in our study will have different hair cut or style, a form of between-subject variability. You can sample a strain of hair of each person or multiple hairs.

1. Fixed effects model: we assume that any systematic variability introduced in our measurement is due to within-subject variability only. In other words, the four men (or women) are more or less equal in terms of hairstyle.
2. Random effects model: this model assumes a more realistic scenario where a man (or woman) is randomly sampled from a population. The model takes into account another factor called "subject", and any unpredictable variation in the data is contributed by an individual or between-subject differences. 

In functional MRI, we employ mixed-effects modeling as accurate statistics when analyzing group data. This is done by taking sample voxel's time series per subject. Theoretically speaking, a fixed-effects model is used when a single subject has multiple runs. In order to assess multiple subjects that have multiple runs, we have to use mixed-effects modeling because every person is understood to be unique and we want to make a statement that can generalize to all people.

Fixed Effects Analysis
Suppose Y_i = X_ib_i + ε_i  is the fixed effects model of our neuroimaging data of N-subjects, each is acquired with n time points. The subject-level GLM consists of p parameter estimates. The within-subject residual ε_i has a multivariate normal distribution with a mean vector 0 and variance-covariance matrix of Σ_i = σ_w².I. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N(n – p). Note how we incorporate the group-level PE from the average of individual subject's PEs, and within-subject error variance from the pooled individual error variance.

Random Effects Analysis
Our group modeling will now have both fixed and random effects. Suppose Y_i = X_ib_i + ε_i  is the model of our neuroimaging data of N-subjects, each is acquired with n time points. The within-subject parameter estimates and residual are b_i and ε_i respectively. In the random-effects model, we further treat the individual subject as coming from a set of the population with its own between-subject variability. Thus, b_i can be further decomposed into b_G + εGi. the group parameter estimates and group error of which the subject belongs. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N – 1. Note that the numerator is the same as the fixed effects model, but the denominator contains more terms. The denominator in random effects is comparatively larger than that in fixed-effects model. So, if you use fixed effects for group analysis, you tend to get unnecessary activation, that is, a false positive.

Notes on Linear Mixed Model
The statistics mentioned in this post are also known by linear mixed model (LMM). With mixed modeling, the model will fit the average intercept and slope as a fixed effect, and each subject is assigned a different intercept/slope. Hence, LMM is an extension of GLM, but this is unlike a standard linear regression model which has only fixed effects. LMM approach is powerful to model data where there is a hierarchical or nested structure, or we deal with multilevel modeling. LMM can also be used in the case of repeated-measures design, or cases when there are some missing data.

Let's take another classic example of LMM: modeling distribution of performance of college students in New York City. Variation in the performance can primarily be due to a random variation at the individual (each student) level. But suppose the city has many colleges. Each college can also be seen as a contributing factor to the performance of each of the individuals at that school, say, due to the good teachers and attitude of the school. Hence those observations cannot be treated as fully independent of each other, but dependent on a higher level grouping called (which college?) - breaking the major assumptions of more traditional linear models. This example reflects different sources of randomness which are in a hierarchy, i.e. individuals-classes-college. This resembles our fMRI dataset, as we move up from an individual subject data to a group-based inference.

LMM is also common for two explanatory variables that are not independent. For example: suppose we want to predict the pitch level as a function of two independent variables: age and gender. Following GLM notation, we can write down:

pitch ~ nationality + gender + error

Notice how age and gender are dependent as coming from the same person or subject. Taking into consideration fixed and random effects with LMM, we rewrite the equation to be:

pitch ~ nationality + gender + (1 | subject) + error

Basic Independent Component Analysis (ICA)

What is ICA?
As an exploratory method, Independent Component Analysis (ICA) provides an alternative to seed-based resting-state fMRI analysis. Traditionally, ICA attempts to solve the cocktail party problem. The story is this. Suppose I have five people with a microphone talking at the same time. How can I separate my speaker output into voices of the first person, second, third, and so on? Indeed, ICA is a powerful method to discover hidden independent patterns or features from a set of data, thus, an exploratory analysis. Unlike GLM, it is model-free because no assumption is made on the shape/pattern of the actual BOLD response.

There are two types of ICA, Temporal and Spatial ICA. The cocktail party problem is an example of Temporal ICA. The data to be separated contains temporal information, the mixing coefficients or weights of which vary uniquely across different locations. In contrast, the dataset in the Spatial ICA contains spatial maps or networks that get activated altogether at a given time. The weights vary uniquely across time for different networks. Spatial ICA is more popular for fMRI data analysis simply because there are more voxels than time points. An image of 64 x 64 x 64 acquired for 250 volumes or TRs contains 262,144 voxels but 250 time points.

Running ICA on the resting-state data will ideally yield to a set of ICs, some of which are clearly related to activation or network, physiological processes (heart rate, breathing), or even some artefacts (e.g. motion, ghosting, slice dropout, noise, etc). Since rs-fMRI is not task-based, we don't have prior knowledge of the temporal waveform of our resting-state networks. This is where data exploration is useful. Spatial ICA attempts to split the data into a set of spatial maps, each with an associated time course. No assumption is made, no experimental paradigm to be specified. The diagram from the FSL website below summarizes this. The content below is thus based on FSL-MELODIC tool.



ICA for Resting-state fMRI
Our observation is the fMRI data. Like standard fMRI analysis, the pipeline starts with the common fMRI data preprocessing steps such as motion and slice timing correction, registration, spatial smoothing, etc. We have to also prepare the preprocessed dataset before ICA. First, simplify the 4D dataset into time vs space. Consider Y_t´m as the whole-brain fMRI image of m voxels and t number of volumes, where t < m. We then mean-center the data by removing the mean spatial map from each row of Y. We normalize the time series variance such that each column of Y will have unit variance, a step called variance normalization. Next, we remove the mean time series from each column of Y. Without variance normalisation, the PCA step below will be biased towards tissue exhibiting high temporal variability.

Fig-1: Temporal std deviation of a sample dataset together with its PCA components (FSL slides).

Two typical yet strong assumptions are made when doing ICA in FSL: (a) the observed source signals are statistically independent; and (b) they are non-Gaussian or not normally distributed. So, ICA is performed on observations that are assumed to be a linear combination of independent sources. Our fMRI image Y is a mixture of independent sources such that Y = A.s , where the A is a t ´ t mixing matrix, and s is a t ´ m spatially independent components.

The first part of ICA is carrying out dimensionality reduction to simplify the problem or computation. This can be achieved with the help of PCA decomposition using SVD. The criterion is that only components that represent a large amount of the dataset remain. The original dataset is now whitened with unit variance and reduced dimension. Rewrite the equation, Y_w = P.Y = P. A.s ; so we have the new mixing matrix A_w = P.A. Note: removing autocorrelation in the fMRI dataset is called whitening, a necessary step in the GLM parametric analyses to make the results more accurate! (Woolrich, et al., 2001). Different statistical software tools (AFNI, FSL, and SPM) have their own steps to accurately model this temporal correlation.

Because Y_w is whitened, its variance-covariance matrix Σ_w is equal to identity matrix I. In other words:  Σ_w = (Y_w Y_w^T) / (n ─ 1) =  (A_ws).(A_ws)^T/ (n ─ 1) = (A_ws s^TA_w^T) / (n ─ 1) = I. This is because the variance-covariance matrix of s.s^T = I as the components are independent to each other. Note that for any orthogonal matrix, A_w^-1 = A_w^T , making the problem solving even simpler.

So now, our unmixing matrix is the inverse of the new orthogonal matrix A_w. The job now is to predict or estimate the inverse of Â_w, such that ŝ = Â_w^-1 Y_w is maximally independent. Note that ŝ and Â_w are estimates since they are both unknown. The criteria of independence can be empirically found by using a certain algorithm, e.g. mutual information, negentropy, etc.

In dealing with fMRI data, specifically, FMRIB Oxford team proposed a probabilistic ICA model in the noisy dataset, a more robust independent component estimation where the noise properties follow is assumed to be Gaussian (Beckmann, et al., 2004). For a more detailed but readable explanation to write this post, I consult Chapter 10 in Ashby, 2011, textbook. I have been using FSL-MELODIC to perform a single-subject ICA and their website is also informative to read.

Theoretically speaking, when all ICs are added together (each one being a 4D signal formed by the outer product of the spatial map and timecourse) they equal the original data. Unlike PCA, ICA enforces independence between the components spatially, while PCA enforces orthogonality both spatially and temporally. Sometimes, the ICs may also share similar time courses. Once all the independent components have been identified, they are ordered according to the degree of importance, that is, % total variance explained.

The last step in the MELODIC is the thresholding stage that produces thresholded ICs overlayed on a background image of your choice. A threshold level of 0.5 (set by default) means that a voxel 'survives' as soon as the probability of being in the 'active' class exceeds the probability of being in the 'background' noise class. This 0.5 assumes we set an equal loss on false-positives and false-negatives. The user guide says that if instead we consider e.g. false-positives as being twice as bad as false-negatives you should change this value to 0.66.

Fig-2: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the physiological signal. The thresholded spatial map illustrates the location of brain regions correlated with the heart rate (~1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.

Fig-3: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the visual and posterior parietal area (0.01 - 0.1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.

By default, FSL MELODIC produces the following files: 

- HTML report with a logfile that collates all steps performed by melodic.

- The mask, the list of ICs, and other statistical information (smoothest : estimated smoothness, PPCA : estimated intrinsic dimensionality estimated from PPCA).

- melodic_mix: an ASCII text file that contains the estimated mixing matrix (in the noise-free case the ICA decomposition is typically written as X = A*C, where X is the original data, A is the mixing matrix, and C is the matrix containing the estimated independent components as its rows). melodic_mix contains #ICs time courses as its columns. Each time course is plotted in the IC report that melodic produces.

- melodic_FTmix: a matrix containing the power spectrum at different freq for the time courses contained in melodic_mix (plotted in the IC report under the time courses)

- Eigenvalues_adjusted : the set of eigenvalues from the initial PCA decomposition (after variance normalisation and adjusting for the dimensionality of X).
 

Group ICA Method
To assess resting-state network with ICA at higher-level or group application, group ICA (gICA) is introduced. One can apply gICA to a set of resting-state data through temporal concatenation gICA; or a set of task-based fMRI using tensor gICA. In temporal concatenation gICA, we are finding common spatial patterns across subjects with unpredictable or different time series. In tensor gICA, we assume that each subject carries a similar time-series pattern, e.g. in task-based fMRI.

Temporal concatenation is recommended for the group-level resting-state pipeline as an alternative method to seed-based analysis. Often, researchers are faced with ambiguity in seed selection and gICA can be useful, especially with the use of dual-regression or back-projection method. For a more detailed explanation, refer to the review paper by V. Calhoun's group (NeuroImage, 2011).

Smith et al., (PNAS, 2009) show that the functional connectivity networks obtained at rest from gICA method are similar to different task-based activation networks obtained from a separate analysis of thousands of subjects. This is quite an important study as the authors show that resting-state networks are not simply random fluctuation, but are recruited when subjects perform the actual task. Several functional networks are identified from resting-state data, e.g. the visual network, sensorimotor network, attention network, limbic network, cerebellar network, the default mode network (DMN), and auditory network. This so-called functional segregation shows that the brain is active even at rest.

Analyzing fMRI data with GLM

What is GLM?
GLM (general linear model in statistics) forms the core of the statistical pipeline in fMRI data analyses. Suppose we want to localize a brain activity associated with moving the left arm. And suppose we utilize a simple blocked design consisting of two alternate blocks of rest and move. We can predict the parts of the brain associated with moving the left arm through a simple prediction. The area in which BOLD signal changes in the same manner as our task-block can be the focus. In essence, what GLM does is multiple regression, i.e. predicting an outcome based on several predictors (regressors, a.k.a independent variables, or explanatory variables). In fMRI, the outcome Y is the voxel time series and this operation is done for each voxel in the brain.

So, what is GLM? There are three basic meanings:

    a)  It is a model to predict or estimate the actual BOLD response based on a known task design.

    b)  It is 'linear' because our prediction is formed through a linear combination of many factors. 

    c)  It is 'general' because we use assumptions to proceed with the statistical tests (t-test, F-test).

Basic statistics remain the same. The variable that we want to predict is the dependent variable. In GLM, the BOLD response of each voxel in the brain is our dependent variable, or to be more precise, it is the voxel time series. The predictors or regressors are also called EVs or explanatory variables. GLM tries to model each voxel's time series as a linear combination of many EVs. 

In other words:

    a)  GLM is a model-fitting strategy done one by one for every voxel, a mass-univariate statistics.

    b)  The contribution of each regressor is called parameter estimate (PE) or β coefficient.
    c)  As in normal regression, the difference between the fitted and actual data is called the residual.


Fundamentally, GLM is a simple linear regression, that is, a model that consists of a predictor (explanatory variable) and a coefficient (slope or β-value). The following diagram sums up what happens when you construct a GLM when you analyze your fMRI data and how it relates to the subsequent calculation of the t-value. A good source for this is the FSL Course slides.

GLM in a Diagram

Fig-1: An illustration of how General Linear Modeling of functional MRI data analysis works.

Key Steps in GLM:

• Our goal is to predict the voxel time series (y) during the task phase using our predicted BOLD response model. Voxel that best predicts the observed time series should be responsible for task performance.
• We assume we roughly know the shape of the BOLD response based on our experimental paradigm. The approach is to employ a canonical HRF or hemodynamic response function such as the double-gamma function. Convolve this with the boxcar function to obtain the model of the BOLD response that corresponds to our blocked design. Note, this function is also called the basis function because you want to express multivariate data using a certain function.
• Determine the whole set of regressors, EVs, or predictors. These include the canonical model mentioned above, motion-related variables, all other noise signals as noise regressors, etc. A collection of EVs is also known as the design matrix. We include noise regressors to take into account artefacts corrupting the actual BOLD response (hence the term: to regress them out).
• A good design matrix should not contain an EV that is highly correlated or a linear combination of another EVs. If this happens, the model is rank-deficient. Our estimates may then fail to represent the actual BOLD response. A well-designed study/task, hence, should ensure that the regressors in the GLM are orthogonal. This topic needs a separate discussion.
• In GLM, the best prediction occurs when the residual is minimized. Thus, to find the best fit of our data is to estimate β values for each voxel such that the residual is minimum.

The solution for unbiased estimates of betas can be found by using the Gauss-Markov Theorem but with some assumptions involved: 

• The residual errors ~ N(0, σ².I),  are normally distributed and independent of each other, 
• The homogeneity of variance is held, 
• Each regressor cannot be a linear combination of other regressors. 

Unfortunately, fMRI observations contain some sort of temporal autocorrelation denoted by a matrix V, making the dirty error variance Σ_e = σ².V. To ensure the data satisfy the assumption, such autocorrelation shall be handled properly. Data preparation such as prewhitening using FILM in FSL or coloring (temporal smoothing in SPM) are meant to achieve this.

Sources of autocorrelation in fMRI dataset include neural source, scanner drifts, physiological signals (respiration and cardiac pulsation), and the head movement. Not accounting for this temporal autocorrelation results in spuriously high fMRI signal at one time point that bleed over to the subsequent time points, which increases the chance of false positive results in task-based studies. The word 'pre' means the operation is added or applied to the original data. The word 'whitening' is a term that is motivated by the fact that white light contains all visible frequencies in equal amounts, i.e. making the spectrum flatter.

What are the steps? 

• First, we fit the GLM ignoring the autocorrelation, estimate the residual errors which contain the correlated voxels and have variance Σ_e. 
• Use the obtained residual to estimate the autocorrelation matrix V. FSL performs a local voxelwise estimation (instead of the global in SPM), which can be modelled quite well with an AR(1) function. Then smoothen this estimate V to correct for any biases using a Tukey Taper, hoping to downweight noisy estimates at higher lags. 
• Once V is obtained, we find a new matrix K such that KVK' = I, where I is the identity matrix. Subsequently, prewhitening is achieved by multiplying both sides of the original GLM equation with this matrix K. 
• We then repeat the GLM again to obtain the correct beta weights. The resulting beta estimates are now unbiased! 

Hypothesis Test: Contrast?

Okay, suppose the whitened GLM is completed. What's next? We want to know which voxels are highly activated due to task performance. How? If the resulting β from the GLM is big and the residual is small, highly likely the voxel is significantly activated during the task. As a quick note: in simple linear regression, β can be treated as a slope. It represents the change in the dependent variable Y resulting from a unit change in the predictor. In statistics, we assess how good the predictors are by comparing whether the PEs are significantly different from zero. Why? If the β = 0, it means the predictor does not contribute anything to the outcome Y.

We are now in the position to perform a t-test on every voxel to know which ones carry significant βs. The formula is given below, the denominator denotes the var[c^T β]. This variance depends on the residual and the critical assumption here is that the residual follows a normal distribution ~ N(0, σ²I).

What is the # dof? On each voxel, we have n time points and n is usually large. A few minutes scan can give you 500 volumes, In our t-test, the degree of freedom = n ─ p. With such a huge # dof, our t-statistics is approximately a z-distribution. Our statistical test gives a map called "Statistical Parametric Map". Because of the number of voxels in the brain, and knowing in reality that each voxel isn't that independent, the problem of multiple comparisons is challenging. Another post will discuss how this can be managed, e.g. through Random Field Theory.

The null hypothesis that a particular voxel isn't significantly related to our task is H₀: c^T β  = 0, where c represents a set of contrast. This contrast is to test the β values and serves as a linear combination of regressors. For p-number of regressors (EVs), the contrast can be written as a p × 1 matrix, each is associated with a beta. In our blocked design, the first regressor is associated with the predicted BOLD response. This allows us to put the first contrast as 1, leaving the rest of unwanted regressors as zeros, c = [1 0 0 ... 0]. This does not look like a normal contrast but as long as there is no singularity in the GLM, the hypothesis test is deemed valid.

Another example, say, if you use two different and independent tasks, e.g. visual and auditory tasks, you want to localize voxels related to your visual, auditory, or both on average. You may then set the contrast to be c₁ = [1 0 0 ... 0], then c₂ = [0 1 0 ... 0], c₃ = [1 1 0 ... 0]. If you want to localize voxels that are significantly more active in the visual than in the auditory task, then c₄ = [1 -1 0 ... 0]. The importance of assigning correct contrasts should not be ignored. Note that we can also perform an F-test, just like ANOVA, to test whether any of the contrast is significant.

Fig-2: An example of a rank-deficient GLM. The design matrix contains two task-related regressors from a blocked experiment (no noise regressors for simplicity). We can equally well use either EV1 or EV2 to explain what we see in the data. In other words, Y can be fit with any linear combination of the EVs, as long as β1 + β2 = 0.9. This yields a problem in hypothesis tests with [1 0] or [0 1] contrasts, but we can still get away with that when testing [1 1] contrast. The computation may still be successful, but your inference hereafter won't be accurate. 

For more detailed explanation on GLM, I find the following sources useful:

      [1]   A nice review article by Martin Monti (Front. Hum. Neurosci., 2011).      
      [2]   Chapter 5, in "Statistical Analysis of fMRI Data", a book by Ashby FG (2011).
      [3]   A more classic Chapter 9, in "Functional MRI: An introduction to methods" by K. Worsley (2001).
      [4]   For this post, I consulted Chapter 7, in "fMRI Techniques and Protocols", Woolrich M, et al. (2009).

Voxel vs Cluster-based Thresholding
In general, the main statistical analysis involves performing GLM at every voxel of the brain. That's why it is also called the mass univariate analysis. Once the statistical test is carried out, we get a z-map or t-map with a large #dof. Each voxel is constructed under the null hypothesis that nothing interesting happens because of our behavioral task. In statistics, the most common (frequentist) approach of a hypothesis testing is to obtain a p-value which is then compared against a certain threshold α. This α also represents the chance of a false positive which is capped at 0.05. Simply put, when we have data with 100,000 voxels, we have 5,000 voxels deemed false positive when we use α = 0.05. The error or bias in drawing a conclusion due to multiple comparison problems is also called familywise error because we are essentially repeating the same t-test for all voxels.

In statistics, the most common correction method for familywise error is Sidak-Bonferroni correction. In fMRI, however, the method is highly conservative. Two more popular methods are widely used in the neuroimaging community: the classic Random Field Theory or GRF (Worsley, K. 1995/1996), the False-discovery Rate or FDR (Benjamini & Hochberg, 1995). A more recent development based on non-parametric statistics using the permutation method is also popular. Detailed discussions about these methods are not presented here.

Post-stats thresholding is the last step in the pipeline to draw the conclusion about our neuroimaging data. This step is done with two objectives in mind. First, we want to know whether a certain voxel or a set of voxels is really active due to the task. Second, we draw that conclusion after having a proper correction for multiple comparisons. Technically, there are two different ways of thresholding:

(1) Voxel-wise thresholding
We can correct for multiple comparisons and do thresholding at every voxel by showing which part of the brain is active at a particular significance level. We reject the null hypothesis that there is no activation if t > u_v. This is called voxel-wise thresholding with threshold u_v. The merit of this method is the high specificity but it faces a serious multiple comparison problem. In earlier days when the scientific community was overly excited about fMRI, the results were reported using voxel-wise thresholding using a more stringent Bonferroni correction, say p < 0.0001, focusing on an area of interest, say frontal motor cortex or visual cortex only. Theoretically, Bonferroni correction is not suitable as it assumes adjacent points are not correlated.

(2) Cluster-based thresholding
Cluster-based thresholding is the more accepted way in the community currently. If a brain region is activated, we expect that not only one, but the surrounding set of voxels gets activated too. Instead of dealing with an individual voxel, we look at a set of voxels called a cluster. Cluster-based thresholding is done in two steps. We first create a binary image (pass/fail) of voxels passing the cluster-forming threshold u_k. Following this, we find another thresholding parameter that is based on either cluster size or peak voxel. For example, we set the cluster parameter with a size greater than a certain value p_k. Here, the sensitivity is typically better but worse specificity. Correction for multiple comparisons is more complicated because it involves a series of contagious voxels. This will be discussed in more detailed in my other post (Random Field Theory, RFT).

Case study: My own data
I'd like to write a simple example of how one can perform statistical analysis (GLM) in FSL. A participant performs a very simple motor task inside the scanner with the right-dominant arm. This is a blocked design with two conditions (rest | move | rest | move |.....) lasting for about 7 minutes (TR=1.69 sec). The scan recorded 250 volumes, meaning, the data would have 250 time points.

The following graphs illustrate what we expect from the GLM analysis to identify voxels that are significantly associated with the task. The red line in the top panel is our observed data, the time series of a particular voxel at a certain coordinate. Ideally, the best fit would see a smooth HRF curve that spans over the 'move' block. The purple plot of a full model fit of the actual time series given by the GLM analysis, that is, given all of the regressors, whether any activation in this voxel can be attributed to any of the task conditions. Lastly, the green line represents only the contrast of interest, and is usually only meaningful when having 3-4 task conditions or taking a simple main effect. Another way of looking at full model vs partial mode fit is this. The former shows what happens when all EVs are used in the fit, whereas the latter shows what happens when only some EVs are used (those involved in the contrast).

The resulting z-map shown at the bottom panel is thresholded at Z = 3.0 with a corrected cluster p-value < 0.05. This simply means that among the voxels that satisfy Z > 3.0, we apply the cluster-based thresholding with a correction for multiple comparisons. In FSL, GRF theory is used as the default correction method.

Fig-3: Outputs of the task-based analysis using GLM. The color map overlayed on the brain image is sometimes called the statistical parameter mapping or SPM (after Friston, et al). Voxelwise temporal autocorrelation must be taken care of by default (prewhitening stage) to ensure that the GLM outputs are valid.

fMRI Data Preprocessing

MRI has become one of the most popular non-invasive brain imaging tools in research. By default, MRI scanners save the data in DICOM format, a standard format for handling, managing, and distributing imaging data. Image data are basically matrices of numbers and can be analyzed further with any software (Matlab for example). We have to do something on the image data before any statistical inference. Why?
- Like many devices, there are sort of noises contaminating the recorded data.
- Participants are not able to stay still ... obviously.
- Our experiment requires them to perform something that may cause the head to move.
- For the hope to increase confidence in the presence of the BOLD signal changes.

There are basically three sets of data one has to acquire when one performs fMRI-related studies:
- The high-resolution anatomical image. This is a T₁-weighted image.
- The functional images in 4D (3D spatial + 1D temporal). These are generally T₂* EPI images.
- Local field information if one wants to take into account field inhomogeneity.

Data Preparation 
The very first step before any processing is to convert DICOM into the format that your software package is using. SPM8 has its own function to import DICOM > NIfTI format. There is one extremely useful conversion software here by Chris Rorden called dcm2nii. Note: People at the MNI prefers MINC format. We can check the post-conversion images by using any MRI viewer software, for example, FSLView. Once you get the files converted, you do some housekeeping. You have to carefully categorize the files into the appropriate folder. For example: T1 anatomical images can be put separately from the multiple runs of functional datasets. 

Once file housekeeping is done, let's preprocess the dataset. To me, the easiest one is to strip the skull from the brain. The skull is not of our interest and it should be removed to simplify the data crunching during the analysis. Apply this step to both the anatomical image and raw EPI images. In addition, you have to remove a few volumes from the raw EPI data. Why? Because the magnetic system requires ~1.5 - 2.0 seconds to reach a steady state.

Head motion detection and correction
Substantial and sudden head motion bring negative effects to image acquisition. First, it contaminates the signal of a particular voxel with the neighboring voxels, making the voxel time series inaccurate. Second, it disrupts the magnetic field homogeneity in the bore that has been adjusted prior to functional scans. Managing head motion should be done before any other signal processing and statistical processes and usually involves the detection and correction or realignment.

Motion detection can be done by first selecting the target or reference volume. The choice of this reference volume can vary. For example: FSL uses the middle volume, AFNI uses either the first or a preselected volume. Often, head motion occurs very abruptly giving rise to spikes or outliers. Such outliers can be identified using a motion outlier script that uses certain metrics, e.g. rms intensity difference between the test and reference volume, or between volume n and n + 1 (DVARS), etc. Once chosen, a threshold to determine an outlier can be found using a boxplot method. Another approach includes a data exploration technique called Independent Component Analysis (ICA) which is more popular for resting-state fMRI. AFNI, though, prefers a different approach called motion scrubbing which means you remove some particular volumes where spikes occur.

Knowing when and how much the movements occur is one thing, but correcting for the contamination is another step. In the "correction" step, all other volumes will go through a 6 DOF rigid body transformation to find the best possible alignment with respect to the reference. The six parameters (or DOF) include 3 XYZ translations and 3 roll/pitch/yaw rotations.

Fig-1: Example of how motion correction analysis estimates the amount of head movement (by FSL-MCFLIRT).

Traditionally, we would want to compute the sum of the squared difference between the target and reference images. We would then minimize this or some cost function through iterative procedures or optimization. Once motion parameters for getting the best alignment have been determined, the new and resampled volumes will be created through spatial interpolation. Why? Because we have to recalculate the new voxel values for each volume after the adjustment. Such "correction" parameters are especially useful as regressors in the next pipeline (statistical analysis using GLM). This is usually a list of 6 parameters according to the rigid body transformation. Often, motion outliers are also used as additional regressors, e.g in FSL.

Slice timing correction
To begin, one 3D image is also called a volume and one scan produces a full functional MRI dataset comprising multiple volumes. When analyzing one 3D image it is assumed that all slices are acquired simultaneously. In reality, however, this isn't the case. Rather, slices are obtained sequentially according to a certain slice order. Thus there is bound to shift in the individual time course across voxels of different slices. For example, imagine we did a scan with a TR = 2.5 sec, i.e. it takes 2.5 sec for a volume to be fully acquired. What this means is that the difference in time between the very first and last slice in a volume would be ~ 2.5 sec. The problem is worsened by the way we acquire the slice, e.g. using an interleaved slice acquisition.

The severity of slice time misalignment depends on the repetition time (TR) and paradigm involved. It is usually very important for an event-related design, which may not be that crucial for blocked-design. All software package that I know of has the correction feature. AFNI and FSL use an almost similar basic method of interpolating the time series. Such interpolation follows the way how the excitation happens in the magnet. It can be ascending or descending or interleaves. FSL provides options to do slice timing correction based on custom-made slice-timing data.

Although controversial, the slice timing correction is usually performed after the head motion correction because the overall effect of the first is less severe than the latter.

Note: Slice time information can be obtained directly from the DICOM header. This can be easily retrieved using a function in Matlab or SPM, but not FSL.

Spatial smoothing (blurring)
This step involves applying a Gaussian kernel to each voxel in a three-dimensional way, which is essentially averaging data points with the neighboring voxels. A Gaussian kernel is a mathematical function that is specified by its width σ (sigma) at half maximum or FWHM, whose center coincides with the voxel center. As a rule of thumb FWHM is selected to be 2x - 3x the voxel size of our functional data, although in practice any value 5-10 mm is very common. Smoothing is also similar to applying a low-pass filter to the original data.

Why do we do spatial smoothing? First, it improves the signal-to-noise ratio. This is related to the assumption that when a voxel is activated, the surrounding voxels are activated as well. Also, we assume the noise of a voxel is not correlated with the noise of the adjacent voxel. Second, it helps to maintain a statistical validity associated with Random Field Theory to solve a multiple comparison problem. Lastly, smoothing is useful during the group-analysis because it improves inter-subject registration by blurring any residual anatomical difference. Risk of overdoing it? The optimal spatial resolution and some desired frequency components are lost if we use too big a kernel. There is also the risk of incorrectly placing the location of the activity or even missing the whole activity itself if it averages out the signal too much.

High-pass filtering
This is basically a voxelwise temporal filter to remove unwanted drift and low-frequency components across time that obscure the actual BOLD changes. This drift is inherent in hardware design and the filtering is common ever since PET is used. High-pass filtering also removes the temporal dc value associated with the image. Common software tool such as SPM uses a discrete cosine function added to the General Linear Model (GLM) to model the drift. FSL uses a high-pass Gaussian filter.

What is the cutoff frequency? FSL gurus recommended that the filter period should be 2x or 3x the task duration per block. A cutoff of 100 sec or 0.01 Hz is therefore very common and also long enough to avoid removing meaningful signals. For an event-related design, there is no clear stimulation period. In order to assess what the cutoff should be, one has to analyze the frequency content of the expected activations.

Some practical examples
Okie, I managed to play around with my own data using different preprocessing steps. The data were obtained from Siemens Trio 3T scanner, with TR=1.69 sec EPI sequence, voxel size 2 × 2 × 2 mm. I stick with FSL5.0 as my software tool. Look how each preprocessing step has an impact on the raw EPI image, both for resting-state and task-based scans. The task is as followed. It is a blocked design between "Move" and "Rest". The subject puts her fist in front of her chest and makes repeated upward movements during "Move", and rests her hand on her chest during "Rest". It's clearly shown that a high-pass filter (HPF) removes low-frequency drift and spatial smoothing with σ = 5.0 mm attenuates the signal amplitude. These outcomes are clearly seen in resting-state data when the brain is not engaged in any special task.

Fig-2: The different outcomes of each preprocessing stage to resting-state data in FSL-FEAT.

Fig-3: Similar treatment but to task-based data, blocked design (Move & Rest, 30-sec each), right arm localizer task. Changes in BOLD signal outweigh the low-frequency drift. In practice, preprocessing is followed by a statistical analysis using the general linear modeling (GLM) framework.

What is the output of any fMRI data analyses? It is the statistical map that shows different activated brain regions associated with the task. What is the effect of different FWHM values on our map? Figure 5 below sums up everything by using FSL software package. Only five selected brain slices are shown through the bottom view:
(a) no preprocessing at all: no slice and motion correction, smoothing, and filtering;
(b) with preprocessing, σ = 0 mm;
(c) with preprocessing, σ = 2.0 mm; also refer to filtered EPI image on the right-hand panel.
(d) with preprocessing, σ = 4.5 mm, but without motion correction (MCFLIRT);
(e) with preprocessing, σ = 4.5 mm; the preferred FWHM is twice the voxel width, not too much!
(f) with preprocessing, σ = 15 mm; also refer to filtered EPI image on the right-hand panel.

Skipping motion correction gives a noisy map because it causes many false-positive clusters. Overblurring or oversmoothing the raw data yields to enlarged activation clusters that are most likely false positive.

Fig-4: Different brain activation maps due to different preprocessing steps. The task is the same as the one in the previous figures. The maps are rendered on the standard MNI 152 template. The M1, premotor, S1, S2, and vermis are shown to be activated during Task execution w.r.t Rest.

Fig-5: Activation maps with respect to rest in coronal, sagittal, and axial views (Z = 3.5, p < 0.05) overlayed on the standard MNI 152 template. Different spatial smoothing parameter is shown with different colors: 2.0 mm (green), 4.5 mm (blue), and 15.0 mm (red).  

  
References
[1]  Some course materials from: http://fsl.fmrib.ox.ac.uk/fslcourse/
[2]  http://support.brainvoyager.com/functional-analysis-preparation/
[3]  Lindquist M. (2008). The statistical analysis of fMRI data. Statistical Science 23: 439–464.
[4]  Ashby, F. Gregory. (2011). Statistical Analysis of fMRI Data, 1st ed. MIT Press.
*** Note: There are abundant piles of past literature dealing with each step of the preprocessing pipeline.