Samir

Investigating the Effects of Smoothing and ICA on High-Resolution fMRI Data

Summary

This project explores how two popular preprocessing methods, smoothing and ICA, influence high resolution fMRI datasets. Using a dataset focused on the brain's motor regions, the motor cortex, the basal ganglia and the cerebellum, I investigated:
1. The effect of smoothing on T-statistics and voxel significance
2. The ability of ICA to delineate task-specific signals from noise

Motivation

The ability to coordinate muscles is fundamental to human motor control, and underlies even simple motions like holding a cup with a hand and flipping a switch with the same arm’s elbow. Past research has not elucidated the brain’s coordination strategy because neurophysiology experiments are confounded by the fact that at a fine spatial scale, motor related neurons correlate with all observable movement parameters [cite] , while at a higher spatial scale, functional neuroimaging experiments have been unable to delineate any generalizable motor organization apart from fractured somatotopic sensorimotor maps [cite].

Probing the motor control regions in a manner that elucidates the structure of the underlying motor controller thus requires novel experiments. High resolution fMRI is ideal for such experiments since it can probe neural correlates of motor activity with a high spatial resolution (a few mm3) and a temporal resolution similar to most motor tasks (a few seconds). The same high resolution fMRI data, however, might be interpreted in multiple ways with minor changes to the processing pipeline. The goal of this project is to understand how two popular preprocessing steps, smoothing and ICA, can influence statistical results.

Methods

A novel motor control dataset was acquired for this study. The data collection started before the course but continued through its duration.

Subject Details

For the data presented in this project, I am including scans for 3 healthy right-handed volunteers (2 males; 19–28 yr of age) from our set of subjects. All subjects were informed about the experiment's details in advance, and gave their informed written consent to a protocol approved by the Institutional Review Board of Stanford University. The subjects were healthy and did not have any psychiatric or neurological disorders at the time or in the past. Subjects were asked to lie supine on the scanner bed, and their heads were surrounded by cloth to comfort and reduce head movements. In addition, they bit down on a bite bar for the duration of the study. The bite bar was customized to each subject's dental structure with a slow-hardening putty. Our controls removed motion related artifacts in almost all our datasets (Fig. 1 shows a representative example).

Figure 1 : The motion estimates and GLM design matrix for a single subject. Note that the mean activity doesn't change with task type and that head motion is not correlated with the tasks. The regressors are orthogonalized and stimulus presentation is randomly interleaved.

Task Details

Subjects were asked to perform writing tasks like drawing a square with a pencil. They were provided visual text stimuli indicating when to plan and execute motions, and when to rest. For instance, one stimulus sequence would be: `Plan : Square' (yellow, 5 sec), `Execute : Square' (green, 8 sec), and `Rest' (red, 3-11 sec, randomized). During the plan phase subjects would plan their motion and possibly adjust grasp position, but would not make any whole arm motions. During the execute phase, subjects would move their entire arm to draw a square, without any finger movement. And during the rest phase, they would place their arm on their torso and rest. Subjects were asked to abandon tasks midway if they could not complete them in time.

The task presentation was randomized and orthogonalized using optseq. A sample design matrix representing one subject's task presentation demonstrates optseq's sequencing (Fig. 1).

Stimulus Presentation

The task description text was displayed on a modified Samsung SyncMaster 305T 30 inch diagonal display (76 cm, 16:10 aspect ratio), built by Resonance Technology (www.mrivideo.com) at a resolution of 1280x800 pixels. Subjects viewed the display through a double mirror at a distance of about 190cm from the head-coil, and the display was inverted so text appeared normal to the subjects. The text was rendered on the screen using VisionEgg, a freely available python based stimulus presentation software.

Scan Sequence Details

We acquired images GE MR 750 3-Tesla scanner and a Nova 32 channel head coil at the Stanford Center for Cognitive and Neurobiological Imaging (CNI). Two series of 510 functional volumes were acquired using a gradient echo, echoplanar sequence with a 134 square matrix, 27 oblique slices, 2mm thick, with a 0.5mm gap between slices. The voxel size was 1.5mm x 1.5mm x 2.0mm, repetition time was 2 sec, echo time was 33ms and flip angle was 75degrees. The slices were adjusted for each subject to include cerebellum, posterior striatum, and motor cortex.

The functional images were overlaid on a co-aligned, high-resolution anatomical scan of the whole brain taken at the end of each session (BRAVO sequence; TR = 7.8 sec; TE = 3.1 ms, flip angle = 12 degrees; matrix, 300x300; 0.8 mm, isotropic).

The Effect of Smoothing

Data

Figure 2 : T-statistic activation maps for subject 1 overlaid on their own anatomical data, for unsmoothed functional (A) and smoothed (3mm isotropic FWHM kernel) functional (B) data. T-values between 3-8 were rendered with a heat map. P<0.001 unc thresholds voxels at t>3. All shown regions also activated at P<0.05 FWE, threshold at t>4.5, in the unsmoothed data. Both functional and anatomical data were coregistered to the mean functional.

Figure 3 : T-statistic activation maps for subject 1 at varying thresholds, for unsmoothed functional (A) and smoothed (3mm isotropic FWHM kernel) functional (B) data. The maps were thresholded with t varying from 0-4 at integer increments and rendered with heat maps. P<0.001 unc thresholds voxels at t>3. P<0.05 FWE thresholds voxels at t>4.5.

Figure 4 : T-statistic activation maps for writing tasks for subject 1 projected on to their own cortical and inflated left hemisphere surfaces. Note how smoothing reduces the number of voxels plotted. T values range from 3-8, and significance is as in Figs. 1 and 2.

Discussion

For the high-resolution data we collected, smoothing substantially decreased the T-statistic across the dataset, and often led to a complete loss of significantly activated voxels in entire brain regions (Fig. 2). While smoothing theoretically increases signal to noise, it did not help our analysis much. The only metric where smoothing held a minor advantage, was that it avoided false positives throughout the data (Fig. 3).

I reconstructed each individual subject's brain surfaces to avoid introducing co-registration and normalization related errors, and created inflated and pial brain models with freesurfer. The projected T-statistics for a subject demonstrate how smoothing dramatically reduced the significance of motor activity (Fig. 4). This effect was conserved across subjects and multiple motion types. Note that projecting on to an inflated brain warps single voxels into circles, making them seem larger than they truly are.

Another important aspect to note is that while anatomical and smoothed functional underlays might display grey matter, CSF related scanner artifacts could have actually created voids at the same position. However, if activation in the void persists after carefully checking the data for motion artifacts and using family-wide error correction you can't reject it. The vasculature is often dense at the boundaries of the cortical surface and if the neural activation in neighboring regions is strongly correlated with the task, then the blood flow in the vasculature might make it seem like the BOLD signal is coming from the artifact region [Personal communications with Prof. Gary Glover, Stanford University].

The task-related variance captured in different components by ICA

Data

Figure 5 : Running MELODIC ICA (part of FSL) on a single session for one subject revealed a complicated eigen spectrum profile, with about 180 predicted significant components. While the first 20 components captured about 40% of the variance across the entire functional data, other components continued to add information.

Figure 6 : While some independent components correlated well with all task regressors, many did not. The components that did not correlate well with any task regressor could potentially be used to denoise the data. Note that all correlations between -0.1 and 0.1 were set to zero to make the plot more readable.

Figure 7 : The p values for whether correlations between independent components and task regressors were significant. Only P<0.005 are shown, others are white.

Figure 8 : The thresholded z-scores for some correlated and anti-correlated components. The correlation analysis helped easily identify motor related and noisy components. The components contained 0.95, 0.86, and 0.77 % of the total explained variance (0.65, 0.59, and 0.53 % of total variance). Note that the cerebellar component actually has an inverted z-score to correlation relationship (+ve z-score = negative correlation).

Discussion

Since smoothing did not seem like a valuable strategy to increase the contrast significance, I next explored whether ICA can delineate noise from task-related activity. The analysis used FSL's standard pipeline, McFLIRT motion correction, brain masking, and 3mm isotropic smoothing. FSL's MELODIC implementation predicted that my data was about 175 dimensional, and produced 175 independent components to describe it (Fig. 5). Manually sorting through the components was not only grueling, but actually counter productive since many components looked like they might be valid "motor components" but did not contain any task related information. Any template based component identification that only looked at spatial localization would suffer a similar fate.

To automate the process of identifying useful components I computed the correlation of each component with all the regressors, essentially applying my GLM model to each component's averaged time series (Fig. 6) and computing the P values (Fig. 7). Doing so revealed the components that contained task-related information. Verifying the components, I found that well correlated components generally appeared next to motor regions (Fig. 8). The fallacy of manually sorting components came out again when the z-score actually involved negative task regressor correlations for a few components. Investigating why MELODIC produced such scores might be interesting for anyone who plans to use its z-scores.

While ICA seems like a promising approach to identify noisy underlying probability distributions in the data, my preliminary investigation doesn't support using it to identify task-related components for high resolution fMRI studies. First, simple GLM analysis give much finer results than ICA's broadly active regions. Second, the thresholded z-scores can not be interpreted on their own and must be compared to some regression model. Third, varying the number of voxels changed the regions in components and their correlations with the task regressors in an arbitrary manner (analyzed data for 7 and 12 components; email me for the results). And finally, the percentage of variance captured by task-correlated components is very small. Moreover, my investigations did not explore the effects of different smoothing levels on the components, which I am sure will be substantial.

Unsmoothed Neural Correlates of Writing Tasks

Data

Figure 5 : T-statistic activation maps for writing tasks for subjects 1-3 projected on to their own cortical and inflated left hemisphere surfaces. T values range from 3-8, and P < 0.05 FWE is at about t=4.5. Notice that while activity is centered around the motor cortex, it is different for all three subjects.

Discussion

With my results, I decided to process my data without smoothing it or denoising it with ICA. Fortunately, since the experiments were well designed and the scan setup at the CNI is really great, I got nice activation correlates which I converted into t-statistic maps using freesurfer. The t-statistic maps for three subjects measuring motor activity during writing displayed significant neural correlates of writing in motor, pre-motor and somatosensory cortices (Fig. 5), as well as cerebellum (not shown). The locations at which voxels correlated significantly with task regressors changed from subject to subject, whose brains had visibly different shapes.

Future Work

These preliminary results are a part of my investigations to ensure that my experimental protocol achieves the highest possible resolution and sensitivity to task-related BOLD signals. In the future, I plan to systematically explore other data-driven algorithms to identify motor features in my data. I also plan to use biomechanical models and robotic control theory to predict the functional organization of the motor regions. Unfortunately, the study is at an early stage and the details will become clearer as it progresses. Please get in touch if you would like to hear more about it.

Appendix

Script

• Code repository with some often used scripts

Software Used

• SPM 8 For the GLM analysis.
• FSL For ICA.
• Freesurfer For creating inflated and pial brains to render SPM's t-maps upon.
• [www.brain.org.au/software/mrtrix/ mrTrix] For plotting SPM t-maps overlaid on the anatomical data.
• Optseq For randomizing the stimulus presentation to ensure orthogonal regressors and minimize subject adaptation.
• VisionEgg For displaying text stimuli to subjects.
• [www.mathworks.com Matlab] For numerous custom processing scripts.

REFERENCES

 [1]