-The Anchorite
==============
In the world of FMRI analysis, certain parts of the brain tend to get a bad rap: Edges of the brain are susceptible to motion artifacts; subcortical regions can display BOLD responses wholly unlike its cortical cousins atop the mantle; and pieces of tissue near sinuses, cavities, and large blood vessels are subject to distortions and artifacts that warp them beyond all recognition. Out of all of the thousands of checkered cubic millimeters squashed together inside your skull, few behave the way we neuroimagers would like.
Two of the cranium's worst offenders, who also incidentally make up the largest share of its real estate, are white matter and cerebrospinal fluid (CSF). Often researchers will treat any signal within these regions as mere noise - that is, theoretically no BOLD signal should be present in these tissue classes, so any activity picked up can be treated as a nuisance regressor.*
Using a term like "nuisance regressor" can be difficult to understand for the FMRI novitiate, similar to hearing other statistical argot bandied about by experienced researchers, such as activity "loading" on a regressor, "mopping up" excess variance, or "faking" a result. When we speak of nuisance regressors, we usually refer to a timecourse, one value per time point for the duration of the run, which represents something we are not interested in, or are trying to divorce from effects that we have good reason to believe actually reflect the underlying neural activity. Because we do not believe they are neurally relevant, we do not convolve them with the hemodynamic response, and insert them into the model as they are.
Head motion, for example, is one of the most widely used nuisance regressors, and are simple to interpret. To create head motion regressors, the motion in x-, y-, and z-directions are recorded at each timepoint and output into text files with one row per timepoint reflecting the amount of movement in a specific direction. Once these are in the model, the model will be compared to each voxel in the brain. Those voxels that show a close correspondence with the amount of motion will be best explained, or load most heavily, on the motion regressors, and mop up some excess variance that is unexplained by the regressors of interest, all of which doesn't really matter anyway because in the end the results are faked because you are lazy, but even after your alterations the paper is rejected by a journal, after which you decide to submit to a lower-tier publication, such as the Mongolian Journal of Irritable Bowel Syndrome.
The purpose of nuisance regressors, therefore, is to account for signal associated with artifacts that we are not interested in, such as head motion or those physiological annoyances that neuroimagers would love to eradicate but are, unfortunately, necessary for the survival of humans, such as breathing. In the case of white matter and CSF, we take the average time course across each tissue class and insert - nay, thrust - these regressors into our model. The result is a model where signal from white matter and CSF loads more onto these nuisance regressors, and helps restrict any effects to grey matter.
==============
To build our nuisance regressor for different tissue classes, we will first need to segment a participant's anatomical image. This can be done a few different ways:
1. AFNI's 3dSeg command;
2. FSL's FAST command;
3. SPM's Segmentation tool; or
4. FreeSurfer's automatic segmentation and parcellation stream.
Of all these, Freesurfer is the most precise, but also the slowest. (Processing times on modern computers of around twenty-four hours per subject are typical.) The other three are slightly less accurate, but also much faster, with AFNI clocking in at about twenty to thirty seconds per subject. AFNI's 3dSeg is what I will use for demonstration purposes, but feel free to use any segmentation tool you wish.
3dSeg is a simple command; it contains several options, and outputs more than just the segmented tissue classes, but the most basic use is probably what you will want:
3dSeg -anat anat+tlrc
Once it completes, a directory called "Segsy" is generated containing a dataset labeled "Classes+tlrc". Overlaying this over the anatomical looks something like this:
Note that each tissue class assigned both a string label and a number. The defaults are:
1 = CSF
2 = Grey Matter
3 = White Matter
To extract any one of these individually, you can use the "equals" operator in 3dcalc:
3dcalc -a Classes+tlrc -expr 'equals(a, 3)' -prefix WM
This would extract only those voxels containing a 3, and dump them into a new dataset, which I have here called WM, and assign those voxels a value of 1; much like making any kind of mask with FMRI data.
Once you have your tissue mask, you need to first resample it to match the grid for whatever dataset you are extracting data from:
3dresample -master errts.0050772+tlrc -inset WM+tlrc -prefix WM_Resampled
You can then extract the average signal at each timepoint in that mask using 3dmaskave:
3dmaskave -quiet -mask WM_Resampled+tlrc errts.0050772 > WM_Timecourse.1D
This resulting nuisance regressor, WM_Timecourse.1D, can then be thrust into your model with the -stim_files option in 3dDeconvolve, which will not convolve it with any basis function.
To check whether the average timecourse looks reasonable, you can feed the output from the 3dmaskave command to 1dplot via the-stdin option:
3dmaskave -quiet -mask WM_Resampled+tlrc errts.0050772 | 1dplot -stdin
All of this, along with a new shirt I purchased at Express Men, is shown in the following video.
*There are a few, such as Yarkoni et al, 2009, who have found results suggesting that BOLD responses can be found within white matter, but for the purposes of this post I will assume that we all think the same and that we all believe that BOLD responses can only occur in grey matter. There, all better.
Thanks for the brief explanation about nuisance regressors. I just entered fMRI world, and the new jargons easily make me confused. You explained it well! Thank you!
ReplyDeleteThis is about six months too late, but you're very welcome! I'm glad it helps.
Delete-Andy
Thanks for this post. I'm wondering if there is any easy way to extract the mean timecourse from WM and CSF with spm8 oder spm12?
ReplyDeletegreetings
dave
Hey Dave,
DeleteThat's a good question. If you've already done the automatic segmentation through the SPM GUI, you should have both white matter and grey matter masks; you could then multiply that with your functional image, and extract the timecourse using code like the following:
V = spm_vol('yourMaskedFunctional')
G = [];
for i = 1:length(V)
G = [G; spm_global(V(i))];
end
figure; plot(G)
I'll have to think more about this, but that is the general idea.
-Andy
Hi Andy,
DeleteJust to follow up on the question asked by Dave regarding using SPM8 or 12 and I'm wondering if you had actually tried or tested it and if you could share more details. Thanks a lot for your help and really appreciate what you are doing!
Jessy
Hey Jessy, I haven't tried that yet; I don't know when I'll get to it in the future, but I'm not as proficient at doing that in SPM as in FSL or AFNI.
DeleteWhat you could do, is try it in AFNI or FSL (using 3dmaskave or fslmeants, respectively) on the masks, and then feed that text file as an unconvolved, "raw" timecourse into the SPM model.
Best,
-Andy
Dear Andy,
ReplyDeleteThank you for this article! I have run the segmentation with SPM and got a probabilistic wm mask. I now would like to create nuisance regressors to remove some scanner drifts that we unfortunately have in the data; I have a feeling that I have to use a conservative wm probability threshold for this, to make sure I do not include any gm/functional activity in the regressor. But I don't know how to select this threshold.. 0.9? 0.8? 0.99? What threshold would you recommend? Or would you recommend do a multiplication, as in your answer to the previous question? (but then some functional activity is still in the mask, right?) thank you!
Hello Irina!
DeleteThe threshold is entirely up to you, and depends on the quality of your data; better image quality will have clearer demarcations between white and grey matter, which will affect your threshold selection. Play around with a few numbers and see which one gives you the most reasonable segmentation of the white matter from the grey matter.
Best,
-Andy
Dear Andy;
ReplyDeleteAll of your tutorial videos are great. Thank you very much. I just started to preprocess rs-fMRI data with FSL and I have learned it through your helpful videos.
As I'm new in rs-fMRI, I am wondering if my understanding is right in regressing out WM and CSF mean signal. I think that I could define ROI for WM or CSF according to known atlas (for instance a sphere with predefined radius as expressed in some articles) then using fslmeants (or MatLab commands), the time series average of ROI can be found and it is subtracted from the data. Would you please kindly help me if I am wrong in this issue?
Thank you very much.
Hi Farzaneh,
DeleteYou could place a sphere, but it would be better to mask out the entire tissue (or liquid) using the segmentation masks. That way, when you average over the entire mask, you get a good general signal for that tissue, which should be different from grey matter.
Best,
-Andy
This comment has been removed by the author.
DeleteDear Andy,
ReplyDeleteUsing FAST and fslmeants, I created 2 text file of CSF and WM signal, but I don't know how I can use GLM GUI to regress out these signals considering that the fMRI data is resting state and no contrast model is used. It would be great if you kindly help FSL users with your great videos.
Best Regards
Hi Farzaneh,
DeleteFirst, select the Custom (1 entry per volume) Basic Shape option on the model setup tab for each type of tissue you are trying to model; then, load your timeseries for each timeseries you have extracted. I think that should do it.
-Andy
Hi Andy!
ReplyDeleteI've been stuck on this for a long time. Its so basic no one ever bothers to explain it. Lets say I have the average WM signal and plug it in to the equation.
Y(BOLD resting data) = X(Mean WM)*B + error
What do I do next?
Hi there,
DeleteIs mean white matter signal the only regressor going into your model? If so, the error time series - i.e., the signal not accounted for by white matter activity - will be what you can then use for resting state analysis. Typically researchers will also include motion regressors and physiological regressors (e.g., heart rate and breathing rate sampled every TR) to account for that variability as well.
Hope this helps!
-Andy
Hah thanks! I just checked back. Still relevant.
DeleteHi, I'm a bit unclear on how to apply the CSF and WM nuisance regressor from this method. I get that I use 3ddeconvolve, but do I run a separate 3ddeconvolve with the two nuisance regressor on the errts file? Or do I use the same 3ddeconvolve cmd I used to get the errts file in the first place, but increase -num_stimts by 2 and of course add in the time coruses with -stim_file?
ReplyDeleteHi Andy,
ReplyDeleteI have segmented T1 image by using FAST and brain extracted by using BET. I have performed preprocessing for motion corrected rsfMRI by using feat. I would like to know how to extract time series of WM and CSF by using FSL?
Thanks
Ramesh
Hi Andy, do you know if people usually normalize the signal in each voxel (i.e., with a z-score across time, or converting to % signal change), before averaging white matter voxels together to make a nuisance regressor?
ReplyDeleteHi Julien,
DeleteYes, the signal should be normalized (usually by converting to percent signal change), so that those voxels are on the same scale as all the other voxels in the brain. The preprocessing streams of all the main software packages (AFNI, SPM, FSL) will take care of this for you, so the residual time series that you extract from should be normalized.
Best,
-Andy
You wrote FIRST for FSL, but you mean FAST for GM/WM/CSF segmentation, right?
ReplyDeleteYes, that's right; thanks for pointing it out. I've corrected it.
Delete-Andy
Dear Andy,
ReplyDeletethanks for this great blog.
I am trying to understand the differences between your suggestions and AFNI's (ex, from afni_proc). I already saw in your answer to another question that you prefer using the average of all data instead of the local voxels in a sphere and I understand it.
Now, what is the difference between creating the mean time course from the errts file as you do, or from the pb03 file as I saw in afni_proc?
(They actually compute the WM stats with pb03, and then they use 3dTprojet on pb04 with the X matrix -- created with 3dDeconvolve while creating the errts file--, using the -dsort option... Getting a little mixed up with all these files!!).
Hey there,
DeleteThe main difference is that 3dTproject is way faster (and can also do bandpass filtering at the same time as it models the nuisance regressors). Other than that, both methods do essentially the same thing.
There's a post on 3dTproject here: http://andysbrainblog.blogspot.com/2014/08/blistering-fast-functional-connectivity.html
Best,
-Andy
Hi Andy,
ReplyDeleteThis is probably a very dumb question. Is it possible to just exclude CSF from analyses, just like how we exclude skull? I am always wondering why ventricles are included in analyses if no signals are expected to stem from there.
Thank you!
Yuqi
Can I first resample the T1 image and then segment the WM and CSF? This will save time.
ReplyDeleteAlso, is the CSF outside the brain region i.e between brain and skull, taken into consideration while regressing it out? Or we just discard that portion of CSF ?
Hi Andy,
ReplyDeleteCould you please elaborate more about how to regress out the noise from the signal? I am doing ICA motion correction. I have identified the motion IC components and have converted their time courses to txt file.
Y=beta0+beta1*IC1+beta2*IC2......+residual
I need this residual to continue doing 1st-level and 2nd-level analysis. But I cannot get the residual images from SPM. I specified the motion IC txt file as "multiple regressors" in "fMRI specification" step and select "Yes" in the following "Model Estimation" for "Writing residual". Then, I got 165 (my scans are 165 time points) separate res_*.nii files. Then, I use them to do 1st-level analysis and SPM keeps erroring.
Many thanks!
Shane