Correcting for Multiple Comparison in fMRI

I like the statements above (source: FSL training slides). In doing statistical inference or hypothesis testing of the fMRI dataset, we face multiple comparison issues. In the earlier post, we say we adopt the mass univariate approach or GLM but we don't know which voxels in the brain the effect is significant. Under the null hypothesis H₀, we say that all voxel activations are merely due to chance. Each voxel has a random value. For the hypothesis testing, e.g. we are repeating the same statistical test with different voxels until the whole brain is covered. If the number of voxels = 100,000 with α = 0.05, then chances are we may encounter 5,000 false positives, i.e. voxels wrongly concluded as being active or related to our manipulation. This basically means that we may get a false positive (Type-I error) in our inference if we don't do a correction properly.

There are some familiwise error (FWE) correction methods. The most stringent and conservative correction of multiple comparisons is the Bonferroni method. There are other variants such as Sidak-Bonferroni, Holm, and Tukey HSD methods. In general, Bonferroni correction is useful if and only if the tests done are statistically independent, that is, it assumes adjacent points are not correlated. This assumption is violated in fMRI data because a strong activation in one voxel is always accompanied by activity in neighboring voxels. Why? In the real brain, the neuronal activation of a point in space is always correlated with the surrounding activation. 

Nichols & Hayasaka (2003) wrote a nice and a bit technical summary of popular methods used to correct for multiple comparisons in the fMRI dataset.

Random Field Theory (RFT)
A nice summary of RFT pipeline is described by Ashby's "Statistical Analysis of fMRI Data" (2011). RFT is applied here to make a statistical inference or hypothesis testing.

(1). What is a Multivariate Normal Distribution?
A univariate distribution is a probability distribution of one random variable. For example: the distribution of height. When the probability follows a normal or Gaussian distribution, it is called univariate normal distribution. If this is applied to a higher dimension containing n-random variables, it becomes a multivariate problem, e.g. a set of variables defining a person such as height, age, arm spread, etc. A multivariate normal distribution is very useful in fMRI statistics because it is used to model the error or residual portion in GLM and PCA/ICA. In a multivariate normal distribution, the relationship among its random variables is a pairwise linear relationship.

(2). What is a Gaussian random field? 
A Gaussian random field (GRF) is an ordered collection of random variables. In GRF, the random variables in the field have a joint-multivariate normal distribution. From our fMRI stats stage using the GLM, we obtain a set of statistical parametric maps, each voxel has its z-statistic. Hence, we are no longer working with the voxelwise image value obtained from the MRI machine, but on the z-scores obtained from the GLM. Under the null hypothesis, these voxels are then treated as a random field. 

Two assumptions are noted while treating the brain voxels as GRF. The first assumption is that the spatial smoothness of the fMRI signal is constant over the brain, and the second assumption is that the spatial autocorrelation function has a specific shape.

(3). What is the basic intuition of RFT?
Suppose we have a 2D slice of z-statistic with 100×100 voxels under the null hypothesis. To fulfill the assumptions, we smooth the data using the Gaussian kernel, a type of spatial filter performed in the usual fMRI preprocessing steps. As a result, instead of the z-map, what we have is resels of 10 x 10. Worsley said that resels are like the smoothed version of the original voxels. They represent the number of independent datasets. What comes next is to determine a cutoff or threshold that will be applied to the resels. Suppose our 2D image which has 100 resels under the null hypothesis is shown below, the Y-axis is the z-score. 

(source: UCL intro slides to RFT).

Let T be a certain threshold applied to the field, such that there are regions that will be above T called the excursion set. The corresponding thresholded image is the Euler Characteristic (χ_T). Intuitively, Euler Characteristic is roughly equal to the # peaks (blobs) with respect to the # valleys (holes, in between peaks) in the image after the thresholding. Importantly, for sufficiently high threshold T we can only see only a few peaks that exist at voxels, where their largest z-value > T. 

Under the null hypothesis, these peaks correspond to # false positives. Hence, the probability of at least observing one false positive is equal to the probability of seeing at least one peak, i.e. χ_T ≥ 1. It turns out that the α_E  is roughly equivalent to the expected value of χ_T , the formula of which has been derived by (Worsley, et al. 1992, 1996). Subsequently, the corrected p-value can be computed through the probability of χ_T ≥ 1 and this is proportional to the volume size and inversely proportional to the smoothness criterion.

(4). Cluster-forming threshold
As with the usual hypothesis testing, we assume that the null hypothesis has to be true for an image voxel to apply GRF to it. Here, we start by forming a so-called cluster step-by-step until we obtain a supervoxel in 3D. We can iterate the process to obtain the distribution of χ_T and z-values. What determines this relationship is the resels only, not the original z-statistic values. If T is the threshold that corresponds to E(χ_T) of 0.05, then by using T, we can expect that any remaining peaks have a probability of 0.05 that they have occurred by chance given that the null hypothesis is true. This α_E value is akin to the Type-I error in the usual hypothesis testing. This Z- threshold value is called the cluster-forming threshold. 

Example: for our 2D resels of 100 resels, the equation gives E(χ_T) = 0.049 for a threshold Z = 3.80, i.e. the probability of getting one or more blobs where z > 3.80 is 0.049. What does it mean? If we threshold our image at Z = 3.80, we can conclude that any blobs that remain have a probability of less than or equal to 0.05 that they have occurred by chance. In other words, this is the case when we have to reject the null. See the illustration above.

For more technical details, this is a good online summary. In SPM and FSL, once we set the cluster-forming threshold (Z = 2.30, with p < 0.05), GRF is automatically applied to obtain the contiguous cluster that survives the threshold while at the same time controlling the Type-I error at 0.05.

(5). How conservative is GRF?
It is good to note that # resel counts is not the same as the number of independent tests. GRF is conservative when the # dof is small, but less conservative than the usual FWE corrections. Without smoothing, the GRF assumes that voxels are independent and the threshold T is equal to the Sidak-Bonferroni threshold, which is a more stringent way of correction. As the kernel size increases, smoothing and spatial correlation increase, resels size increases, and so T decreases. The choice of the most optimal smoothness places a serious problem in applying GRF to the fMRI dataset. Another thing to note that in a 3D brain volume, the resel takes into account the shape of the volume. A more sophisticated method using search volume can be used (Worsley, et al. 1996).

False Discovery Rate
FDR approach is totally different from controlling familywise error. Rather than limiting the errors V themselves to occur, we limit the proportion of significant results that are false positive, that is, V/R. Thus, FDR is the probability of V/R and the corrected p-value guarantees that FDR < q. The goal of FDR is finally to find regions with p-values associated with q smaller than 0.05, say.

What is the algorithm in finding q?
1. Convert N number of z-values to p-values. Let the value of q = 0.05.
2. In the case of activation, p- is the area under the null distribution to the right of z
3. In the case of no activation, p- is the area under the null distribution to the left of z.
4. Rank or order all p-values from the smallest to the largest. Let the k-th smallest p-value as P[k].
5. The z-value of the k-th smallest p-value is significant if P[k] < qk/N;  otherwise it isn't significant.

Compared to Bonferroni and GRF thresholding, FDR almost always indicates more and larger significant regions in the brain. But this is also because FDR does not control for Type-I error α_E, but rather, the rate of its occurrence as a whole.

Permutation-based Method
The permutation-based method is considered non-parametric as it doesn't depend on the nature of the distribution and correlation that exists in our image data. Recent opinions support the permutation-based approach as an alternative to GLM-based statistics. To see why this makes sense, the task-based block design will be used. A participant performs a task in one block and rest in the other. If a specific location of the brain is required for task performance, we should expect a task-related activation during task block as reflected in the BOLD response. However, this activation does not appear during the rest block.

Now, we reshuffle or scramble the data across time points or TRs. If a particular region is involved in the task, we should not expect any strong correlation between the predicted BOLD response and the task performance because the relationship has collapsed. The "reshuffled data" behaves just like a random noise. Reshuffling does not change the mean and variance of the dataset, that is, we don't tamper the data. Although attractive, the method has one drawback: the computation time is much longer as you have to keep randomizing the time points. Although reshuffling the data for 500 or 1000 times is quite a norm, there is no consensus of how many times this should be repeated.

How can this method correct for multiple comparison problem? It sounds simple, but I think I have to read more. Therefore, a separate blog entry will attempt to discuss this method with my own stuff to play. Anyways... what it does in correcting for multiple comparison is mentioned below:
1. Randomly shuffle the data point (TRs). Reorder every voxel time series in the ROI so that the
    TRs agree with the shuffle.
2. Run the GLM to produce a statistical parametric map. Find the maximum z-value (z_max) of the map.
3. Repeat Step 1 & 2 so that eventually we have M maps and thus z_max. Usually M = 500 or 1000.
4. Rank the order of maximum z-values, the largest being z_max[1] down to z_max[M].
5. The corrected threshold T = z_max.(Mα_E + 1), where α_E is the familywise error.

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.

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.

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.

How should I do this?

I cannot catch up! My work requires me to master stuffs from three different domains: neuroscience, motor behavior, and statistics. It is simply too tiring to catch up with the amount of information that I have to learn with the time spent in writing this down online. From now on, I will just write a point of interest that I recently came across. It will be updated and revisited again and again in the future if I have more findings.

This is especially true when, for example, someone asked me about a few questions that I'm not too sure, but there is no sufficient time for me to find all of the answers. Furthermore, it is better not to retype the summary or basic theories of a certain topic. Citing the appropriate source here directly is more practical and time-efficient!

Okie, the 3rd year of my PhD officially has started :(

From Neuroanatomy to Cognition

White Matter Fibres
The white matter was briefly mentioned in an earlier post, so this is sort of a continuation of the brain's gross anatomy. The white matter is located underneath the cortical gray matter and composed of fatty myelinated axons. It is an integral part of the central nervous system that transmits messages very rapidly. It basically has 3 types of fiber bundles: the projection fibers, commissural fibers, and association fibers.

1. Projection fibers are bi-directional, afferent, and efferent bundles. They appear as radiating bundles in the white matter that exit the cerebral cortex and converge towards the brainstem. One bundle carries visual information through the optic radiation. Near the subcortical nuclei, these axons form a compact band known as the internal capsule with anterior and posterior limbs. Afferent (sensory) fibers: mainly the thalamocortical bundles going to the various region of the cerebral cortex. The efferent fibers of the internal capsule arise from the cerebral cortex. They form various tracts, e.g. corticothalamic, corticobulbar, corticospinal, and corticopontine bundles. 
2. The axons part of the corpus callosum forms the commissural fibers. At different callosal segment, they have different connections: the rostrum (orbitofrontal), genu (frontal lobe), body (sensorimotor and posterior parietal), and splenium (posterior temporal and occipital). Other commissural fibers are the anterior commissure, connecting the olfactory system bilaterally.
3. The association fibers form the bi-directional cortico-cortical bridges connecting areas within the same hemisphere. They can be classified as short and long fasciculus:
        - Superior longitudinal fasciculus connects frontal and parietal lobes.
        - Occipito-frontal fasciculus connects frontal and occipital lobes.
        - Arcuate fasciculus connects the frontal with posterior temporal lobes.
        - Uncinate fasciculus connects orbitofrontal with anterior temporal lobes.
        - Inferior longitudinal fasciculus connects temporal and occipital lobes.
        - Extreme capsule fasciculus connects lateral temporal and lateral frontal lobes.

Of interest is the coronal section of the cerebral hemisphere from the insula moving inwards to the thalamus. The external capsule connects the motor cortex to the putamen and is unidirectional. The internal capsule connects specific thalamic nuclei to the specific cortical area and hence it is bidirectional.

The most common way to study the white matter is through MRI which can be observed well on T1-weighted, T2-weighted, and FLAIR sequences. More recently, scientists become more interested in modeling brain development over puberty and brain decline associated with aging. Fun facts: Gray matter volume increases in early childhood but declines after puberty. However, white matter volume progressively increases over time, supporting the concept of neural plasticity.

Cerebral organization
The cerebral cortex is organized into six layers that arise from the time of its development. This is the characteristic of the neocortex. Only the piriform cortex and the hippocampal formation, the oldest cortical structures phylogenetically or paleocortex or allocortex, do not exhibit this six-layer arrangement. The projection fibers are more deep-rooted, while the association and commissural fibers are more superficial. Three principal types of cells found in the cortex include the pyramidal, stellate, and fusiform neurons. Their fibers are arranged either tangentially or radially across layers.

Pyramidal cells, with a shape of a triangle with the top end going up to the surface (apical) and the horizontally running dendrites (basal), constitute the most in various cortical layers. The axons are either going down to the white matter (as projection fibers) or to other cortical areas (as association fibers). The biggest pyramidal cell, the Betz cell, is found only in Layer V of the precentral gyrus or motor cortex. Unlike pyramidal cells, granule or stellate cells are small, polygonal or triangular in shape. They are found in all layers, but especially numerous in Layer IV. Fusiform neurons are spindle-like cells found mostly in the deepest cortical layer, their long axis going vertically upward. Apart from these three types of cells, we encounter others, e.g. horizontal cells found mostly in the superficial layers. The works of Cajal and Golgi are crucial in deepening our understanding on these cells.

Fig-1: Six different cortical layers of the cerebral cortex, layer-I being the most superficial.

In brief, six-layered architecture can be described as follow:
a). Layer I (molecular layer), has few cell bodies, mostly axons, Layer II (external granular layer).
b). Layer III (external pyramidal layer), cells forming mainly association or commissural fibers.
c). Layer IV (internal granular layer), mainly the incoming afferent fibers from the thalamus.
d). Layer V (internal pyramidal layer), mainly efferent projection fibers.
e). Layer VI (multiform, fusiform layer).

What is the relationship between this architecture with the earlier functional lobes? Layer III plays a major role in cortico-cortical connections. Layer IV is predominant in sensory areas in the parietal and temporal lobes, e.g. the postcentral gyrus. These regions are granular. Layer V, on the other hand, is predominant in motor areas, e.g. precentral gyrus.

Fig-2: The distribution of different cortical composition: (1) Agranular; (2) Granular - frontal (dysgranular); (3) Granular - parietal; (4) Granular - occipital; and (5) Koniocortex. Only cortical motor areas are agranular.

Principal neurotransmitters

A variety of neurotransmitters is associated with neurons of the cerebral cortex. Among those, we have glutamate, aspartate, and γ-aminobutyric acid (GABA). Pyramidal cells are the main efferent neurons that are predominantly glutaminergic and are excitatory. Most interneurons within the cortex, however, are GABAergic and are inhibitory. They are bridging the afferent and efferent fibers together. Therefore the outputs of the cortex are modulated by a variety of cortical afferents via interneurons. 

A variety of neuropeptides or monoamines are also found in the cerebral cortex; they influence not only populations of neurons but also local metabolic activity and vascular smooth muscle. The most important monoamines in the cortex are (1) norepinephrine, which originates from the locus ceruleus of the pons and distributes sparsely to all cortical layers; (2) dopamine, which arises from the substantia nigra–pars compacta and the adjacent ventral tegmental area and is found in moderate amounts in layers I and VI and sparsely in layers II to V; and (3) serotonin, which arises from the raphe nuclei and distributes heavily to all cortical layers.

Cognition and the Brain
The study of human cognition and the brain is the heart of a classic science popularly known as neuropsychology. The interests existed since the time of Descartes, Gall, Broca, and so on, who studied the link between a neurological condition (e.g. lesions) and certain behavioral or psychological processes. A classic theory, phrenology, says that the brain is divided into discrete and unique areas responsible for a particular function only. The mastery of certain skills can be deduced by the bigger skeletal landmark of the head. An opposing view at that time held that there is no localization of brain functions and that the functions (what they called "Mind") are distributed across different parts of the brain. With more discoveries, modern neuroscience later thought that the brain is divided into many functional specialization. For example, one may use fMRI to elucidate brain areas associated with some behavioural tasks. One fundamental characteristic of the central nervous system is parallelism, that is, a large number of functions are simultaneously processed along two or more pathways. As a result, the damage of one pathway can allow other pathway to function, and that one brain function can be performed not only strictly by one area. 

Modern neuropsychology enjoys a multidisciplinary collaboration among cognitive scientists, physiologists, neuroscientists, and clinical psychologists. Originally, the field drew strong attention when Paul Broca came into contact with a patient undergoing a progressive speech disorder in 1861, who could only produce "tan". After the patient died, Broca found out that his inferior frontal gyrus (IFG) was damaged. Named after Broca, the type of such behavioral deficit linked to the damage of IFG is then called Broca's aphasia. Note: IFG is rostral to the mouth/orofacial musculature of the cortical motor area.

Fig-3: The difference between Broca's and Wernicke's aphasia together with affected areas on the left hemisphere.

Broca's finding was further developed with the findings of Carl Wernicke. He found that in a certain type of language disorder, the patients were able to produce speech but unable to comprehend the conversation. Called Wernicke's aphasia, the damage is found to be around the posterior part of the superior temporal gyrus (STG). This aphasia is not equal to deafness, for the person with Wernicke's aphasia is able to detect sound but unable to make sense of it. He further hypothesized that there is a link between IFG and STG and this is crucial in language. To be able to converse well, one has to first listen and understand the sentences one hears. Note: STG is near to the primary and secondary auditory cortex.

In the 1870s, John Hughlings Jackson proposed that the cerebral cortex is organized hierarchically and that some cortical areas are for higher-order functions (or cognitive) that are neither fully sensory nor motor. These brain areas are called association areas because they serve to associate sensory inputs to motor response and conduct mental processes related to sensorimotor behavior. The mental processes that Jackson attributed to these areas include interpretation of sensory information, the association of perceptions with previous experience, focusing of attention, and exploration of the environment. Jackson's finding is supported by clinical works. The major helps come from surgical rooms of patients with damage or lesion on the specific are, or people with underlying conditions. Other methods include experimental studies with monkeys and rats and the use of non-invasive brain imaging technology.

Before ending, I wish to mention major associative areas in the human brain important in cognition:

1. The posterior association area: the margin of the parietal, temporal, and occipital lobes. It integrates information from several sensory modalities such as vision, space, and body senses. It is also involved in language. Separate studies by Holmes and Luria on wounded soldiers found that bilateral injuries to the posterolateral parietal lobe yield to normal visual acuity but the soldiers were unable to scan visually or reach for an object of interest. They could not process together with the visual information when asked to describe in words what that they saw. This shows that the region is critical for integrating different sensory modalities and for using that integrated information to direct behavior. 
2. The anterior association area: the prefrontal region, rostral to postcentral gyrus. It is involved in the planning of action, shaping behavior, and judgment; a more popular term is the "Executive function". The most popular case showing how the injured prefrontal region leads to behavioral problems is perhaps of Phineas Gage. A series of clinical tests, e.g. the Tower of London test and the Wisconsin Card Sorting Test (WCST), can be used to diagnose people with neuropsychological disorders who have lost their executive functions, such as schizophrenia. WCST is primarily considered a test of executive functions, particularly abstract reasoning and cognitive flexibility in response to external changes.
3. The limbic association area: along the lower medial end of the cerebral hemisphere. It is for emotion, learning, and memory. Its involvement in learning and memory comes from the well-known study on patient H.M. by B. Milner in 1960s after both medial temporal lobes had been removed. She first demonstrated the remarkably selective role of this part of the brain in converting short-term into long-term memory. Studies in monkeys have helped establish that association areas in the medial temporal lobe, including the hippocampal formation, receive information from virtually every other association area. In other words, the hippocampal formation is able to sample the whole stream of ongoing cognitive activity and thereby relate different aspects of a single event so that they can be recalled as a coherent experience.

More recently, cognitive neuroscience is recognized as another separate field, combining neuroscience, neurophysiology, and psychology. Scientists now agree that the three areas (the triad) of executive function are working memory, flexible thinking, and inhibitory control.

Sensorimotor Control & Learning (Part I)

"Motor learning", the main theme of my current lab, lies at the intersection between motor behavior and neuropsychology of learning. It helps us to understand how a person acquires and learns new movements, e.g. to dance, play golf, or adopt a new language. The theme is studied depending on the organ where the voluntary movements are produced: the arm, legs, jaw, and eyes.

The current summary is based on modern studies from the 1990s by a group of engineers and computational scientists, outside the domain of neuropsychology and kinesiology. Emphasis will be made on the upper limb (arm reaching) as its core manipulation.

Components of Sensorimotor Control
First, gathering sensory information associated with the task. When one wants to move or perform a task with one's limb, visual information is gathered through saccades. This gaze behavior is also task-specific. Our brain appears to be able to filter irrelevant sensory inputs. For example: the study in inattentional blindness when one fails to notice prominent visual stimuli unrelated to the task one is attending. It is worth noting that sensory streams are temporally delayed and noisy.

Second, motor tasks involve a sequence of decision-making processes in the presence of delay and noise. Why does the noise come into play? Because both sensory and motor system are inherently noisy, arising naturally at the molecular, synaptic, and system levels. The final product is that trial-by-trial movement production to the same target in space is bound to exhibit some variability. These underlying risks can be viewed with the context of reward, especially during a period of motor learning. Even when a person faces the same motor learning task, one can be a risk-averse (exploitation) or risk-seeking (exploration).

Third, our nervous system can be modeled as a controller. Traditionally, there are two basic controllers for motor control and learning: the feedback and feedforward controller. Feedback control, as the name implies, refers to the control of voluntary movements using sensory feedback. In contrast, feedforward control doesn't depend on any feedback mechanisms. Given that our sensory inflow has an inherent delay (150 - 200 msec), it becomes unreliable to rely on for an accurate movement control. To meet the demand of the motor tasks, we often rely on feedforward control which is a predictive control. As we produce certain movement, we make also make prediction about the sensory consequences of that movement. This is done using of efference copy of the motor command and the difference between the predicted sensory consequence and the actual sensory feedback will be used in state estimation.

The last controller is related to the biomechanical properties of the body and the tools used. It is within the topic of arm impedance. Impedance control depends on a few factors such as arm stiffness, i.e. how springy the musculatures are. Like the internal model, impedance control is also inspired by concepts in engineering and biomechanics. Example: we can modify the way we grip a tool (hand stiffness, arm stiffness) produced through co-contraction of the opposing muscles. Although co-contraction can be a solution to the motor task, it is inherently unstable as the sensorimotor system is noisy.

Recent Techniques or Methodologies
We can study motor learning in various ways, It can be studied through behavioral studies and quantitative movement analysis. Modern motor learning literature typically includes three well-known behavioral paradigms:

1. Sequence learning: a type of motor skill learning where it employs serial reaction time tasks (SRTT). Here, a participant is asked to make a sequence of button pressing, key tap, or finger flexion. Motor performance is measured by the reaction or response time. It does not directly deal with the kinematic and dynamic features of motor learning. In more specific ways, learning is measured by the difference in reaction time between the random sequence and learned sequence.
2. Visuomotor rotation: this can be achieved by e.g., providing a set of prism worn by the participant or by a certain mechanism to distort the association between the visual feedback and the actual arm movement. The performance is measured by the movement deviation. This method introduces a mismatch between two related sensory inputs: visual and proprioception.
3. Force field paradigm: a participant performs reaching movement with a robotic manipulandum, capable of producing a velocity-dependent force that perturbs the movement trajectory. The presence of the force changes the dynamic of the motor task. The process of reaching motor performance signifies adaptation. The sudden removal of the force causes the trajectory to deflect to the opposite direction know as an after-effect. A more novel idea probes trial-by-trial performance in terms of the magnitude of the lateral force the participant produces. This is achieved by introducing catch trials in the form of force channels. This is the method used by Shadmehr and colleagues.

Lesion studies and non-invasive stimulation (TMS and tDCS) are able to complement the methods. Recently, neuroimaging methods are employed to learn regions of the brain associated with the behavioral tasks involved. Scientists employ engineering and computational modeling to represent the brain as a system or controller. In error-based learning, for example, motor adaptation can be captured by a linear time-invariant model (LTI). With this framework, in each trial, we learn new movements by employing an optimization algorithm (e.g. Kalman Filter).

Types of Motor Learning Processes
The processes of motor learning can be classified by the type of information the motor system learns.

1. Error-based learning: motor learning through the presence of an error, i.e. the discrepancy between the desired trajectory and the actual movement outcome. The term error also means the mismatch between the predicted sensory consequences and the observed sensory feedback. Three important behavioral paradigms to study this type of learning include visuomotor rotation, prism goggle, and force-field adaptation. Error reduction happens reasonably quick and this type of learning is known as adaptation. Improvements in adaptation reach plateau after 8-10 trials. We will focus more on this type of motor learning processes as it has been widely studied for the past decade.
2. Reinforcement learning: normally observed in a redundant system. This type of learning can happen even when there is no error involved. Learning is achieved through exploration by finding the best solution in the solution manifold. It is highly dependent on the available rewards, e.g. points, punishment, currency. A study by Izawa & Shadmehr (2011) shows that reinforcement learning, in some circumstances, can substitute for adaptation when there is uncertainty about, or no, sensory prediction error.
3. Use-dependent learning: not strictly a learning process, but rather, adaptation through repetitive movements. Movements to a certain direction are able to reduce variability in that direction and induce a bias towards this trained direction when reaching to other directions.
4. Learning by observation: typically it involves watching others doing the movements. This type stems from the findings of mirror neurons. Observational learning may include learning from predicting error by observing the action of others.
5. Structural learning: learning to extract common features of different task variants. When we know the underlying structure of the task, learning can be faster, e.g. learning to swing a tennis racket bears similarity with learning using a squash racket.

What is Internal Model?
Modern research of human motor behavior has been marked by the incorporation of control engineering theories, in particular, the concepts of the internal model (see: Jordan, 1995; Kawato et al. 1987). A controller G_c(s) is used to control the process G_p(s). A good controller should be able to represent the process to be controlled. It is said that the motor system is composed of the limbs (i.e. the plant) and the controller in the nervous system, the internal model.

The internal model is an approximation of the inverse dynamics of the system being controlled. It is a model that mimics the behavior of the natural process being controlled, which refers to our motor system. This model can be adapted at any time to a novel environment, making it a suitable computational model for motor learning (refer to studies by Shadmehr's group). There are 2 variants of the internal model: forward model and inverse model.

1. Forward model predicts sensory consequences from the efference copy generated during movement. Forward model is likened to a motor-to-sensory mapping. The efference copy is issued in conjunction with the motor command from the CNS. The model tries to anticipate the next state so that the movement goal is achieved and the error is minimized. 
2. Inverse model tries to approximate motor commands through an inverse transformation from the incoming sensory streams. This model is reactive rather than predictive. There is a close relationship between model (1) and (2).

The internal model is used to explain motor adaptation. As mentioned, it involves a decrease in sensory prediction error through trial-by-trial adjustments in the forward model. Accordingly, the update of the forward model is translated into an update of motor commands. Mathematically, the internal model is well captured by linear time-invariant (LTI) state-space models, which have sensory errors or perturbations as inputs, sensorimotor mappings as hidden variables, and the learned or adapted motor commands as the output.

Why can't we tickle ourselves? When we tickle our body, the central nervous system predicts the sensory consequence using the efference copy of "tickling". At the same time, there is this somatic sensation generated by the "tickling". As this sensation matches the predicted sensory consequence through the forward model, the comparator circuit in the CNS doesn't detect any mismatch.

Can Motor Learning Generalize?
After going through training of a task in one context or situation, a person is able to perform as well to a similar task but in a different context or situation. This concept is called generalization. When generalization is beneficial, it is usually termed transfer. Conversely, when it is detrimental, it is termed interference. Traditionally, the studies of generalization made use of dynamic or force field paradigm. The principles derived from those studies are associated with the concept of the internal model. On top of that, generalization is related to another concept called "motor memory". If interference occurs, the transfer of learning fails.

Using force field paradigm, it is thought that motor adaptation is able to generalize in the intrinsic coordinate system, i.e. based on internal muscular patterns of activity. A salient example of intrinsic transfer is when writing "9" by right and left hand. On the other hand, using the visuomotor paradigm (Krakauer et al., 2000), motor learning generalizes in the extrinsic coordinate system, i.e. based on the external spatial coordinate frame. Transfer in motor learning has also been studied in relevant to transfer across different movement direction. There is a limited transfer of dynamic (Gandolfo et al., 1996; Sainburg et al., 1999).

Transfer occurs in different configurations of the same arm (Shadmehr & Mussa-Ivaldi, 1994; Ghez et al., 2000; Malfait et al., 2002; Shadmehr & Moussavi, 2000). How about the interlimb transfer? Tranfer occurs from the dominant arm to the non-dominant arm and this happens in the extrinsic coordinate system (Criscimagna-Hemminger et al., 2003). The opposite is not true. Further, the interlimb transfer from the dominant to non-dominant hand occurs only when the force field is introduced abruptly (Malfait & Ostry, 2004). The gradual force field, on the other hand, does not cause the apparent interlimb transfer. It seems that interlimb transfer is regarded as a cognitive process.

The Concept of Motor Memory
The initial part of learning involves more cognitive processes, where one makes use of one's memory buffer to carry out and finish the task. This readily available, temporary buffer or space is known as working memory. A popular example of working memory is when you solve mathematical problems. The later part of learning is the period when motor performance stabilizes and involves consolidation, a term related to long-term storage of motor skills. This is why after a year of not playing the piano (or skiing), we are still able to play it as well. However, what is stored inside the memory (e.g. motor commands, task dynamic, somatic experience, etc.) is still debatable and the nature of consolidation is also conflicting, e.g. Caithness et al (2006).

Memories that can be consciously recalled are named declarative memories, e.g. memory of events, words, or facts. Conversely, memories on skills and knowledge to perform some particular actions are called procedural memories, e.g. the ability to walk or ski. Typically, the domain of motor learning deals with procedural memory, the nature of which is interesting to characterize. Smith et al. introduced a two-rate state-space model of the force-field adaptation. The model says that adaptation consists of fast learning with poor retention and slow learning with more stable retention. Krakauer & Shadmehr discuss whether the formation of such memory progresses over time from a labile state, which is susceptible to interference to a stable state, which is resistant to such interference.

The memory from experiences obtained from a rapidly changing environment leads to faster decaying and unstable motor memory. As opposed, more stable memory is achieved when exposed to a gradually changing environment. What happens when, after learning A, a person immediately learns B? Called retrograde interference, task B is able to disrupt the consolidation of A. This phenomenon does not appear after a longer period of training. There is even evidence suggesting sleeps enhance/improve consolidation. Sometimes, although we forget to do a certain task, a quick relearning is sufficient to meet the expected performance. Such a phenomenon is called saving, that is, faster relearning.

Lastly, the concepts of implicit and explicit processes have a place in the context of motor learning. Explicit processes require declarative knowledge of something. When a person learns to make a golf swing, voluntary explicit processes include, for example, adjustment to the weight of the golf stick or the knowledge on the target location. In contrast, proficiency in skill performance itself or the correct timing (or speed) involves implicit processes. In some conditions, implicit planning may override explicit strategies during a visuomotor adaptation task (Mazzoni & Krakauer, 2006). Scientists are still debating which of the two are dominant during motor learning, at different stages of learning.

The famous case of HM who couldn't recall practicing a mirror writing task but performed well in the task several days later prompts scientists to think that motor learning is purely implicit. Adaptation is seen as an implicit process, but it does not rule out the involvement of explicit processes. Such "cognitive" explicit processes are thought to serve as a form of learning strategy. Taylor and Ivry (2011) found two competing processes: explicit knowledge of target error and implicit knowledge of sensory prediction error (a la the usual adaptation mechanism). Keisler and Shadmehr used an interesting approach to examine declarative memory contribution to force-field adaptation. Subjects were adapted to force-A and then a brief exposure to force-B. After a 3-min interval, they experienced channel trials. At the same time, they had to memorize words in between. This memorization interfered with the memory of the second task B.

References  
[1]  Wolpert D.M., Diedrichsen J & Flanagan J.R. (2011). "Principles of sensorimotor learning". Nature Rev. Neurosci. 12: 739-751.
[2]  Krakauer, J. W. and P. Mazzoni (2011). "Human sensorimotor learning: adaptation, skill, and beyond." Curr  Opin Neurobiol, 21(4): 636-644.