Journals: How to write or read?

As a new grad student, one may wonder where to start if you either (1) write a new manuscript; or (2) read a journal in a lab meeting or school assignment. There is no hardcoded way to answer this. My supervisor uses the following steps in writing a fresh manuscript:

(1) Write the method section: the most straightforward way. You know what you did the most!

(2) Continue by writing the result section.

(3) Continue to bring the discussion, highlight the meaning of the main findings.

(4) Write an introduction and abstract the last.

Writing the introduction and abstract is preferentially the last because having gone through the previous three steps, one is able to have a much better overview or picture of the "value" one wishes to contribute. Note that I don't use the word "meaning" or "scientific contribution" because even if a paper has a good contribution but fails to convey its value to a broader audience, it fails to be recognized.

Reading a published journal is another thing to learn. In our lab meeting, we usually start with the abstract to get a rough idea the paper intention. Then we move on method section and figures (yeah pictures and diagrams, graphs, etc.) We then read the discussion in greater details. Nobody cares much about the introduction.

Sensorimotor Control & Learning (Part II)

Limits of Adaptation Paradigms
In laboratory settings, motor learning is often studied in the context of motor adaptation paradigms, in which subjects must learn to compensate for a systematic perturbation, either in the through manipulating visual feedback (Krakauer et al. 2000) or a change in the dynamics of the arm manipulandum (Shadmehr and Mussa-Ivaldi 1994). In the first case, visual feedback is altered resulting in a change in movement kinematics. In the latter, change in dynamics is created through the introduction of the opposing force field which is perpendicular to the movement direction. In motor adaptation, what is typically observed is a monotonic improvement in performance that is initially rapid, and then slows to an asymptotic level close to initial baseline levels of performance. The progress of learning is well described by an exponential fit, implying that the amount of improvement on each trial is proportional to the error (Thoroughman and Shadmehr 2000; Donchin et al. 2003). According to these adaptation paradigms, motor learning is predominantly mediated by a specific mechanism that is based on changing an internal forward model.

In real life, we may see the motor adaptation in the case, for example, when one learns to adapt when holding a heavy tennis racket or counteracting the fluid dynamic under-water. However, not all motor learning falls under the definition of motor adaptation. For example, when we learn to synthesize entirely novel movements in the absence of perturbation, adaptation fails to explain these learning processes. Performance improvements take shape from being incompetent as a naive learner to full proficiency. However, such improvements are far slower than in adaptation paradigms: while tens of trials are usually enough to reach an asymptotic level after perturbation, performance in these more complex tasks continues to improve over hundreds of trials or even across a few days.

Haith & Krakauer (2013) define this long-term reduction in movement variability as skill learning and argue that such learning is associated with incrementally improving the quality of one's movements with practice. Skill learning has been studied in the laboratory setting using a few tasks, e.g. maneuvering a cursor along a constrained path (Shmuelof et al. 2012) or through a series of via points (Reis et al. 2009). overall variability in task performance reduces substantially over days of practice, even though subjects immediately exhibit near-perfect performance at slow speeds.

Learning a "Skill"
Skill is another dimension of motor learning. Krakauer et al. defined skill as a shift of speed-accuracy trade-off function when there is no systematic perturbation to movement trajectory. This can be measured before and after a series of training protocols, e.g through tracing an arc or reaching a point in space. Skill learning involves a slower process, or improvement and it is distinct from adaptation. Adaptation to error is not a skill because the performance reaches plateau once the trajectory returns to the baseline trajectory, i.e. there is no further improvement possible. Krakauer mentioned that performance improvement can occur through:
     (a) Better state estimation (improved forward models, or processing of sensory feedback);
     (b) Better motor execution (improved signal-to-noise ratio in motor output).
Which one is the limiting factor is unknown.

In light of skill learning, an obvious question is this: how can we model skill learning? As mentioned above, error-based motor adaptation has some limitations in real-life cases. Some studies have posited that adaptation also cannot fully explain motor learning processes such as saving, and learning in the absence of error signal. In Huang et al. (2011), subjects experienced visuomotor perturbation. Following adaptation, they went through a sufficient washout block to reset the adaptation. After imposing the perturbation the second time, there is no indication of faster relearning (saving) than the initial learning of naive subjects. However, with the introduction of reward feedback during a task with the same paradigm in another group of subjects, saving is observed following a sufficient washout block.

Model-based and Model-free?
Much of the work on error-based learning uses an internal model to explain motor learning. This is also called model-based learning because it requires building, updating, or adapting an internal model through performing the task. It turns out that another process is possible to support motor learning, called model-free learning. In this model-free type, the goal is to learn directly through a process of trial and error, to explore the space of potential actions in each state, and keep track of which states and actions lead to successful outcomes (rewarded). Indeed, reward-based tasks (reinforcement learning) is an example of model-free learning. Huang et al (2011) and Krakauer (2013) posited both model-free and model-based adaptation are working hand-in-hand to cause saving or faster relearning.

Building an internal model as learning progresses is a salient feature of the model-based approach. Once found, learning can be generalized to multiple but similar task structures. This, however, bears heavy computational costs. In contrast, we usually talk about control policy in reinforcement learning, or model-free learning. A control policy here means the selection of a single action per trial, or describes an ongoing stream of motor commands in continuous time according to the instantaneous state.  No intermediate model representation and no calculation required to transform a forward model into motor commands in the control policy. Model-free learning tends to deliver superior performance on a particular task in the long run because they do not rely so heavily on noisy computations each and every time a movement must be made. But the disadvantage is that the learning scope is rather restricted to the task performed during training. Even if the reward structure of the task changes in a known way, one must start from scratch (or at least from some previous but incorrect control policy). This is in sharp contrast to the flexibility offered by model-based learning.

Optimal Feedback Control
While in daily lives, the internal models require constant adjustment especially in the context of learning or adaptation. Todorov and Jordan proposed that the key to such adjustment is the presence of feedback. There are primarily two types: own sensory feedback (proprioception, vision) and the predicted sensory consequences from the internal model (efference copy).

The computation framework that links internal model, feedback, motor cost and reward, and optimization is called the optimal feedback control (OFC) (Todorov and Jordan, 2002) championed by scientists such as Todorov, Kording, etc. It has been regarded as a comprehensive theory of motor coordination in redundant systems, such as our joints that have multiple degree-of-freedom. The theory says that, for a given motor task there is this optimal control policy to update the motor plan using a cost function of effort and accuracy. An OFC uses an optimal estimate of the state of the system, generated through sensory feedback and efferent copy, and uses this feedback to adjust its output towards a specific goal.

In OFC, a cost function of effort and accuracy terms can be optimized, assuming that unbiased estimates of some key parameters (e.g. for forward dynamic model) are available a priori, to derive a feedback control policy for a given task goal. It is thought that during motor learning, a person first learns to adapt to the task dynamic, not the task-irrelevant goal. This step is then followed by the optimization of a cost function to do the task. In other words, there is no clear one-to-one mapping between tasks and actions in the brain. Instead, the naive learner will aim to "converge" into the best solution.

How does OFC relate to skill learning? What is supposedly occurring in the training period is a set of process that leads to better performance: convergence on the optimal policy, or improved execution of the control policy itself, perhaps through an increased signal-to-noise ratio via expanded neural representations. Either of these possibilities could be the explanation for shifts in the speed-accuracy tradeoff, and reductions in variability described in motor skill learning studies. Ultimately, the authors suggest that the process of converging can be both model-based and model-free.

Lastly there are several theories suggesting the neural substrates underpinning motor control and learning according to this framework:

1. Basal ganglia (striatum) helps to monitor the reward and cost of the motor commands generated (reward-based learning).
2. Cerebellum helps to predicting sensory consequences and monitor mismatch, otherwise known as error monitoring.
3. Parietal cortex combines the expected sensory consequences and the actual sensory feedback, a process analogous to  state estimation, a place for multisensory integration and the body awareness.
4. Premotor and primary motor cortex assign the feedback gain for the visual and proprioceptive states, transforming the belief about the states into the actual motor commands downstream.

References  
[1]  Haith and Krakauer (2013). "Model-based and Model-free Mechanisms in Human Motor Learning". Adv Exp Med Biol. 782:1–21..
[2]  Krakauer and Mazzoni (2011). "Human sensorimotor learning: adaptation, skill, and beyond." Curr  Opin Neurobiol, 21(4): 636-644. 

[3]  Shadmehr & Krakauer (2008). "A computational neuroanatomy for motor control". Experimental Brain Res, 185: 358-381. 

Notes on Resting-state fMRI Analyses (Part II)

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

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

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

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

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

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

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

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

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

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

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

 

Slice Timing Correction?

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

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

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

Adding additional EVs to GLM

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

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

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

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

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

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

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

Denoising nuisance components with ICA

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

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

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

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

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

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

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

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

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

Notes on Resting-state fMRI Analyses (Part I)

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

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

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

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

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

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

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

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

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

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

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

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