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.

Functional Connectivity - a brief overview

Background
During its initial application in scientific research, fMRI was used primarily with task-based behavioral studies. It was used to locate brain activation associated with a certain task or stimulus. For example, we want to study V1 in response to seeing patterns of a checkerboard, or S1 in response to a tactile stimulus on our left forearm. Since Biswal et al., (1995) found a temporal correlation in the sensorimotor area in the absence of tasks, the interest in the so-called spontaneous fluctuations in the brain has accelerated. This type of fMRI, the resting-state fMRI, is unique because the frequency content is low < 0.1 Hz. That is why it was initially thought of as mere physiological noise.

In another study, Fox et al. (2005) reveal a set of interesting brain regions where they get deactivated after a person engages in a task. These regions are known to routinely exhibit task-negative responses (deactivations) during attention-demanding tasks. For terminology sake, brain regions that get deactivated during task performance are called task-positive. When talking about spontaneous fluctuations, we have the terms correlation and anticorrelation regions. Anticorrelation signifies a form of temporal negative correlation. In that study, the fluctuation that is negatively correlated with those in the attention and executive systems is later called the default mode network.

Statistical maps produced by this spontaneous BOLD acquisition is thought to reflect functional connectivity network or resting-state network (RSN). Some authors (Seekey, 2007) used the word intrinsic connectivity to denote two distinct and dissociable networks, the "salience network," (dACC and orbital frontoinsular cortices, with connectivity to subcortical and limbic structures), and "executive-control network" that links dorsolateral frontal (dPFC) and parietal neocortices of the brain. These intrinsic networks are correlated with task-based MRI measured.

Fig-1: (A) BOLD activity during finger tapping and spontaneous fluctuation with a seed located at point a. (B) Spontaneous deactivation observed during task performance using PET. Regions that are negatively correlated with the seed a is called the default mode network.

Methodological Perspective
As mentioned, the only difference between rs-fMRI and task-based fMRI is that the acquisition is performed at rest. The subjects typically lie down in the scanner with their eyes open and looking at a cross-hair, or with eyes closed. Some differences in brain activations exist between the two variants. The minimal requirements of an rs-fMRI study make it easy to scan a wide variety of populations, subjects may get bored and fall asleep. A recent study indicates that falling asleep is the major challenge for subjects and that fixating subjects rarely fall asleep.

The way we analyze resting-state fMRI data is different from the method used in the task-based fMRI. This is because we do not have any clear task model for our design matrix. Having said that, some steps in the preprocessing pipeline are still valid here. For example, the removal of non-brain tissues and skulls, spatial smoothing, high pass filtering, image registration, and motion correction.

Resting-state fMRI suffers two essential drawbacks. First, separating wanted signals from the image artifacts produced by magnetic field distortion (e.g. near the sinuses), signal dropout, and motion can be challenging. Unless we can correct or model these events, the presence of artifacts potentially reduces the temporal correlation between two brain areas. Image distortion due to sinuses can be corrected using field map correction. FSL/FreeSurfer software package has a subroutine called BBR that helps to correct this when registering the EPI to T1-image.

Second, physiological signals such as cardiac and respiratory cycles can also modulate the wanted signals via several mechanisms (Murphy et al., 2013). More recent works have shown that, with a very high sampling rate or TR < 400 msec, we can segregate the frequency contents of spontaneous brain signals. But such MRI sequences are still work-in-progress. There are, however, a few methods to correct for this:

• Regress out the average time series obtained within the white matter region and ventricles, which is mostly of no interest.
• Regress out respiratory and heart-rate signals acquired from the actual physiological recording.
• Regress out global (whole-brain) average time series. This method is still debatable as it allegedly introduces the systematic anticorrelation areas.
• Use ICA to identify artifactual components around the perimeter, cardiac-related, or anything with the frequency > 0.1 Hz.
• Use band-pass filtering to limit the frequency content to be 0.01 - 0.1 Hz.

Once we have more or less 'clean up' our dataset, there are two most common and essential methods in resting-state data analysis, i.e. seed-based and ICA (Independent Component Analysis) methods. The methodological aspects of rs-fMRI are massive and it is always evolving. But some key points I'd like to write about RSN analyses:

• The frequency spectrum is usually between 0.01 - 0.11 Hz, i.e. the 'low frequency' spontaneous fluctuation, although a wider range up to 0.25 Hz was recently accepted.
• There are two most commonly used methods to determine RSN: the seed-based correlation analysis (SCA) and Independent Component Analysis (ICA).
• The SCA method is a model-driven analysis that requires a strong apriori region of interest selection. This method allows us to use the usual framework of general linear modeling (GLM). Example: if we are studying the effect of playing tennis on the brain, we should select M1 as one of our seeds.
• ICA is a data-driven analysis and does not require apriori assumption. ICA is also very attractive in identifying various neural activities at rest, collectively reflected as independent components, e.g. physiological artifacts and jerky head movements.
• It is important that changes in RSN do not constitute the causality of our behavioral paradigm.

Exploring Brain Organization
In several studies, ICA has been used as a robust method to segregate different functional networks of the brain (RSNs). With improvements in the MRI sequences, more and more details can be picked up using rs-fMRI techniques. The number of wanted components was found to be 10 of 25 in (Damoiseaux et al., 2006), 45 of 70 in (Smith et al., 2009), and ~23 of ~200 in (Marcus et al., 2013). Another method is more like a seed-based approach, which involves calculating pairwise correlations between all voxels and using clustering algorithms to identify groups of highly correlated voxels. Those studies identified the brain as being partitioned into somewhere around 15–20 large-scale clusters (e.g. Power et al., 2011; Yeo et al., 2011). Refer to Fig-2.

Fig-2: Various resting state networks revealing different functional areas of the brain (motor, sensory, visual etc).

Some of these networks have been known anatomically for long ago but others, such as attentional networks, are new. Interestingly, in some cases, functional connectivity grouped brain areas that were not so widely recognized as connected. With sophisticated data-driven approaches, research has been spent to refine boundaries between adjacent functional brain regions.

... Resting-state correlations grouped a set of regions and made it easier to recognize that they shared a variety of specific characteristics, bolstering the case that these regions form a functional system. The large-scale (system-level) patterns in resting-state activity, therefore serve as a useful organizing framework for interpreting results and patterns in other modalities.

Recently, a new paradigm of analyzing rs-fMRI data called the network-based approach emerges. When viewing the data as a network, less focus is given on the properties of a single brain area, but more within the larger neural system. However, it is wise to note that a network in this sense is not based on anatomical or physical neuronal networks. Instead, it is from time-varying signals arising of either BOLD (functional networks) or DTI (diffusion), which are essentially mathematical models with some underlying assumptions, that may not fully reflect the physical entity.


Source: Power, JD, Schlaggar BL, Petersen, SE (2014). "Studying Brain Organization via Spontaneous fMRI Signal", Neuron-Primer. Some other notable papers: (Raichle, 2010) for a historical and metabolic perspective; (Deco et al., 2011; Hutchison et al., 2013) for dynamical perspectives; (Bullmore and Sporns, 2012; Sporns, 2014) for network perspectives; (Murphy et al., 2013) for a methods perspective; (Lee et al., 2013) for a clinical perspective; and (Buckner et al., 2013; Craddock et al., 2013) for general perspectives.

Functional MRI - a brief overview

How does BOLD-fMRI work?
To begin, two material properties are important in BOLD, i.e. diamagnetism and paramagnetism. In essence, a diamagnetic material does not introduce a significant change in the magnetic field, whereas a paramagnetic material tends to increase the magnetic field. If these two types of material are close to each other, they cause a local distortion of the magnetic field near the interface. The field becomes less homogeneous. Brain tissue is mainly diamagnetic. In contrast, the magnetic property of the blood may change depending on the oxygen molecules attached to haemoglobin. This is crucial. When the blood contains more haemoglobin without oxygen attached (deoxyhaemoglobin or deoxyHB), it is paramagnetic.

The more deoxyHB the blood has, the more local field distortion it creates. The local field inhomogeneity causes faster spin dephasing in the transverse plane, causing a lower T₂* value. In this case, the image intensity drops. Now, what happens when there is neuronal activity? More oxygen molecules are needed by the neurons or brain tissues, so the oxyHB concentration increase and deoxyHB concentration drops. With more oxyHB the blood becomes less paramagnetic, causing the image intensity to increase. These properties are being exploited in fMRI to capture neuronal activities in the brain. The EPI sequence has been known to be highly sensitive to such changes in magnetic properties, making it the most popular fMRI method to use.

The way the BOLD behaves in the event of neuronal activity is called the haemodynamic response, or BOLD response. The physiology of this response is not straightforward and depends on, e.g. the cerebral blood flow (CBF), the cerebral blood volume (CBV), and the metabolic rate of oxygen consumption (CMRO2). In response to a stimulus, the CBF goes up to deliver more oxygen to the site of neuronal activation. On the other hand, the CMRO2 is increased or more oxygen is consumed, which reduces the BOLD effect.

Fig-1: The relationship between a stimulus, neuronal activity, neurovascular coupling, and BOLD in fMRI scanning [2].

 

Does fMRI measures brain activity? No. It does not directly measure neuronal activity, but rather, it uses blood deoxyHB level as a proxy or indirect measure of neuronal or functional activation. In other words, fMRI measures the degree of neurovascular coupling. Scientists have noted that while the neuronal activation is very fast, the BOLD response is slower.

Experimental Paradigm using fMRI
The experimental paradigm is directly related to research questions in mind and influenced by the fact that the BOLD response is slow. Although the response is more or less reproducible, the shape and the onset may vary depending on the brain region and stimulus duration. Refer to the diagram above. A good paradigm is able to take into account the slow response but carries high statistical power for making any conclusion. At the same time, it should also ensure that the task is not biased, and prevents subjects from anticipating or getting bored, that is.

1. Blocked design: by far the most common paradigm in functional MRI. In this case, one block represents one task or experimental condition, and one scan session involves more than 1 block. The duration may range from 20 - 35 seconds, allowing a fully restored or complete profile of the haemodynamic response (HRF). The HRF can be viewed as a filter (Josephs & Henson, 1999). The most efficient design is a sinusoidal modulation of neural activity with T = 25 sec (e.g., boxcar with 12 sec on/ 12sec off), capturing fully the BOLD signal and its peak. We should design the block in that way. The signal of one particular block is then compared with the haemodynamic signal produced during a rest or baseline period. Thus, the blocked design is actually a subtraction or a contrast between Task vs. Rest brain activity. We can always expand this by using more tasks within a scan (e.g. other stimuli or conditions) which we would compare against the REST block. If any, the interaction effect between task conditions must be taken into account. With regards to this kind of design,
• Advantage: simpler in execution, high statistical power, does not require an accurate HRF model.
• Disadvantage: doesn't allow separation of individual trials, induce boredom and anticipation, and is not suitable for all behavioural tasks.
2. Event-related design: this is the second paradigm where individual events related to the different tasks or experimental conditions are measured. Here, an event is presented at a certain short duration with inter-trial stimulus (ISI) time, and is assumed to evoke a set of neural responses in the brain. The task presentation does not follow a block-by-block arrangement but is presented in a random fashion, each may last only for 2-3 seconds. Event-related design requires the MRI pulse sequence to be fast enough to catch up with the changing task event (e.g. with a relatively shorter TR), giving a higher temporal resolution. The advantage of this paradigm is its flexibility in the experimental design, and the tendency to prevent boredom or fatigue. In practice, there are a few variants such as rapid ER design, jittered ER, and randomized ER. 
• Advantage: flexible, remove anticipation, can separate response to different stages.
• Disadvantage: tedious implementation, low statistical power and sensitivity, require good HRF model (sensitive to error), thus requiring more #trials per stimulus.

Fig-2: The difference between blocked and event-related design with three different behavior conditions.

An important finding that makes the event-related paradigm simpler is the fact that the BOLD response of the event tends to be evoked similarly even when the response of the event before that has not decayed fully. In other words, the responses sum up linearly.

fMRI Signal and Noise
In fMRI, the evoked BOLD response is our signal of interest whose behaviour is not straightforward. Scientists have spent efforts to model this response because this is the first step before making any inferences. It allows us to know in the time domain which one is activation, which one is not. The most common model is the one that assumes the BOLD response to be a linear time-invariant system. Under this assumption, there is a linear relationship between neuronal response to a stimulus and the BOLD response. It is time-invariant and does not depend on any previous stimuli. With this assumption, its characterization is known in a noisy system.

Using this framework, a burst of neural activity or spike can be presented as discrete impulse responses. Then, the observed BOLD response of a voxel can be modelled as the convolution between the incoming stimulus waveform and the impulse response. The resulting response is now called the canonical haemodynamic response function or simply, HRF. The general agreement is to use the double gamma function as the HRF. What are the drawbacks of this model? The linearity assumption may be too simplistic. Also, the shape and onset of BOLD responses may vary across subjects. Notably, the same region doing different functions for the same task may show different evoked responses. Scientists have proposed more robust models for HRF (see [1] and Glover et al., 1999).

The fMRI signals are prone to corruption due to noise and artifacts, collectively known as nuisance signals. Just imagine! The signal change is usually about 2% of the total signal magnitude. Unwanted signals can typically be in the form of:

• Hardware noise: thermal noise (higher magnetic field strength gives more noise) and the scanner drift (usually f < 0.01 Hz, we can filter this out or model it).
• Participant's head movements. Sometimes it appears as a sudden spike.
• Physiological noise: heart rate and respiration, the most challenging one to model/remove. The spectral components of heart-related noise are between 0.9 - 1.0 Hz, while respiration, 0.3 - 0.4 Hz.
• Others: structural-related noise. In 2007, Fox et al reported that spontaneous BOLD follow a 1/f distribution (pink noise), meaning that there is increasing power in the low frequencies.

The presence of nuisance signals distort the wanted BOLD response. This eventually obscures the actual neural activations seen in the image. In other words, noise reduces detection sensitivity.

Brain Connectivity and fMRI
fMRI is extremely useful to identify brain areas associated with a particular task. But what happens if we scan the brain at rest? The brain is never at rest. In 1995, Bishwal found that there is spontaneous low-frequency fluctuation of BOLD in the human brain at rest, that is when the brain is not engaged in doing any specific tasks. The term "resting-state" became popular. Separate research by Raichle and colleagues found specific brain regions called the default mode network or DMN. The unique feature of this network is that the activity decreases when the subjects are engaged in tasks. A group of scientists from FMRIB-Oxford, has identified several consistent RSNs such as those of the visual cortex, sensorimotor, and executive function.

There are a few reasons why this resting-state fMRI (rs-fMRI) is attractive. It does not require the subjects to perform any task inside the scanner. This is important if the devices are not MRI compatible. Until recently, there is an increasing number of publications showing the application of rs-fMRI in the clinical setting such as Alzheimer's Disease, ageing brain, epilepsy, and some pharmacological studies.

In short, fMRI is useful to localize brain activities and study the brain at rest. Recently, scientists become more interested in studying how more than one location interact in the brain. Terms such as network and connectivity are then introduced. We now have three different types of connectivity:

• Anatomical connectivity: as the name implies, it is a hardwired structure of one brain region with the other. It can be studied elegantly with Diffusion Tensor Imaging. Such a technique complements earlier histological methods such as retrograde tracing.
• Functional connectivity: connectivity of two or more brain regions whose time series are correlated. This is studied using resting-state fMRI performed while the person is at rest. Applying ICA on resting state data produces different sets of spatial maps called the functional connectivity network.
• Effective connectivity: connectivity of two or more brain regions where one region influences the other regions. Rather than depicting temporal correlation, effective connectivity emphasizes causal relationship among brain areas.

Data Analysis Pipeline
I would like to end this post by presenting the most common data analysis pipeline. Whatever innovative pipeline one takes, it bears the same objective: to allow valid statistical inferences. There are two categories of data analysis pipeline: task-based and resting-state data. The most common statistical framework used in the analysis is called GLM, the general linear modelling. Here, we input a certain design matrix or schema associated with the experimental or task paradigm and find the brain regions that fit the schema the most. A certain post-hoc (correction) step is required to keep the statistical principles valid. The model-free method called the Independent Component Analysis (ICA) is more attractive to work with resting-state images.

Fig-3: The most common data analysis pipeline in functional MRI. The preprocessing steps are more or less fixed, but the researcher has to choose between the model-based (GLM) or model-free (data-driven) method. Task-based analysis mostly employs GLM, while resting-state fMRI employs a data-driven ICA method. Although there are many versions to this, the basic idea remains the same.

 

Once the experiments have been conducted, we obtain a series of imaging data in DICOM format. The data have to go through preprocessing steps before we perform any statistical analysis. The pipeline produces a final product as a statistical parametric map (after Friston), a form of graphical representation where we can visually see parts of the brain associated with our experimental paradigm. NOTE: There is no one correct pipeline that fits-for-all scenarios. Most researchers tailor it to their research needs.

References
[1]  Jezzard, P., Mathhews, P. M., and Smith, S. M. (2001). Functional MRI: An Introduction to Methods. Oxford University Press.
[2]  Arthurs, O.J. and Boniface S. (2002). How well do we understand the neural origins of the fMRI BOLD signal?. TRENDS 
      Neurosci, vol. 25: 27-31
[3]  Lindquist M. (2008). The statistical analysis of fMRI data. Statistical Science 23: 439–464.
[4]  Cole, D.M., Smith, S.M., Beckmann, C.F. (2010). Advances and Pitfalls in the Analysis and Interpretation of Resting-State FMRI
      Data. Front Syst. Neurosci., vol. 4.