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.

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.