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.