Introduction to Beta Series Analysis in AFNI

So far we have covered functional connectivity analysis with resting-state datasets. These analyses focus on the overall timecourse of BOLD activity for a person who is at rest, with the assumption that random fluctuations in the BOLD response that are correlated implies some sort of communication happening between them. Whether or not this assumption is true, and whether or not flying lobsters from Mars descend to Earth every night and play Scrabble on top of the Chrysler Building, is still a matter of intense debate.

However, the majority of studies are not resting-state experiments, but instead have participants perform a wide variety of interesting tasks, such as estimating how much a luxury car costs while surrounded by blonde supermodels with cleavages so ample you could lose a backhoe in there. If the participant is close enough to the actual price without going over, then he -

No, wait! Sorry, I was describing The Price is Right. The actual things participants do in FMRI experiments are much more bizarre, involving tasks such as listening to whales mate with each other*, or learning associations between receiving painful electrical shocks and pictures of different-colored shades of earwax. And, worst of all, they aren't even surrounded by buxom supermodels while inside the scanner**.

In any case, these studies can attempt to ask the same questions raised by resting-state experiments: Whether there is any sort of functional correlation between different voxels within different conditions. However, traditional analyses which average beta estimates across all trials in a condition cannot answer this, since any variance in that condition is lost after averaging all of the individual betas together.

Beta series analysis (Rissman, Gazzaley, & D'Esposito, 2004), on the other hand, is interested in the trial-by-trial variability for each condition, under the assumption that voxels which show a similar pattern of individual trial estimates over time are interacting with each other. This is the same concept that we used for correlating timecourses of BOLD activity while a participant was at rest; all we are doing now is applying it to statistical maps where the timecourse is instead a concatenated series of betas.

Figure 1 from Rissman et al (2004). The caption is mostly self-explanatory, but note that for the beta-series correlation, we are looking at the amplitude estimates, or peaks of each waveform for each condition in each trial. Therefore, we would starting stringing together betas for each condition, and the resulting timecourse for, say, the Cue condition would be similar to drawing a line between the peaks of the first waveform in each trial (the greenish-looking ones).

The first step to do this is to put each individual trial into your model; which, mercifully, is easy to do with AFNI. Instead of using the -stim_times option that one normally uses, instead use -stim_times_IM, which will generate a beta for each individual trial for that condition. A similar process can be done in SPM and FSL, but as far as I know, each trial has to be coded and entered separately, which can take a long time; there are ways to code around this, but they are more complicated.



Assuming that you have run your 3dDeconvolve script with the -stim_times_IM option, however, you should now have each individual beta for that condition output into your statistics dataset. The last preparation step is to extract them with a backhoe, or - if you have somehow lost yours - with a tool such as 3dbucket, which can easily extract the necessary beta weights (Here I am focusing on beta weights for trials where participants made a button press with their left hand; modify this to reflect which condition you are interested in):

3dbucket -prefix Left_Betas stats.s204+tlrc'[15..69(2)]'

As a reminder, the quotations and brackets mean to do a sub-brik selection; the ellipses mean to take those sub-briks between the boundaries specified; and the 2 in parentheses means to extract every other beta weight, since these statistics are interleaved with T-statistics, which we will want to avoid.



Tomorrow we will finish up how to do this for a single subject. (It's not too late to turn back!)



*Fensterwhacker, D. T. (2011). A Whale of a Night: An Investigation into the Neural Correlates of Listening to Whales Mating. Journal of Mongolian Neuropsychiatry, 2, 113-120.

**Unless you participate in one of my studies, of course.

Extracting and Regressing Out Signal in White Matter and CSF

A man's at odds to know his mind because his mind is aught he has to know it with. He can know his heart, but he don't want to. Rightly so. Best not to look in there.


-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 FIRST 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.

AFNI Command of the Week: cdf

Not necessarily a neuroimaging-specific tool, cdf simply converts between p-values and t-statistics (or F-statistics) using the cumulative distribution function. Supply the test that you did, followed by the t-statistic (or p-value) and degrees of freedom, e.g.:

cdf -t2p fitt 3.4 15
p = 0.00396 #A t-statstic of 3.4 with 15 degrees of freedom yields a p-value of 0.00396

cdf -p2t fitt 0.001 30
t = 3.65 #We would need a t-statistic of 3.65 or greater to reach a p-value of 0.001

Degrees of freedom can be found by using 3dinfo on a statistical dataset, and then looking at the value of "statpar" associated with your statistical test of interest. Degrees of freedom can also be calculated as the number of time points minus the number of regressors in your model; for example, if you have 1200 time points and 40 regressors, then the degrees of freedom will be 1200-40 = 1160. In the following X-matrix (generated by the command "aiv X.jpg"), the first 25 columns represent regressors accounting for any drift during that run; the next nine columns are the regressors of interest; and the last six columns are motion regressors.

 

The degrees of freedom in neuroimaging data can be a little tricky to interpret, as FMRI time series are temporally autocorrelated; in other words, the value of one time point can be predicted, to a degree, by neighboring timepoints. Therefore, using ordinary GLM estimation techniques can lead to an inflated degrees of freedom. To rectify this, instead of using 3dDeconvolve, use 3dREMLfit, which will attempt to account for this autocorrelation.


AFNI Command of the Week: 3dZcutup and 3dZcat

When I used to work at OSU, on the lab wiki I would put up a new AFNI command every week, detailing a program that isn't necessarily used all that often, but has some interesting applications for the user looking for more ways to manipulate their data. I plan to do the same on this blog, in the hopes that someone might find them useful.

One such tool that came to my attention a couple of weeks ago was 3dZcutup, a program for taking apart individual slice or groups of slices, in order to rearrange them or, more commonly, to perform statistical analyses on only one slice at a time, if computer memory becomes an issue. The usage is simple: Supply an input dataset, a prefix for your output dataset, and specify the range of slices you want to dump into the output dataset. For example, say you have a functional dataset r01+orig with 35 slices in the z-direction; if you wish to output only the first half of the slices into one dataset and the second half of the slices into another dataset, you could do something like the following:

3dZcutup -prefix bottomHalf -keep 0 16 r01+orig
3dZcutup -prefix topHalf -keep 17 34 r01+orig


Recall that the slices start at slice 0, which is why the last slice in this dataset is labeled 34. The output datasets for these commands would look something like this:

TopHalf

BottomHalf


In order to rearrange these slices, either to recreate the original dataset or to inverse the slices, you can collate the slices with the complement to 3dZcutup, 3dZcat:

3dZcat -prefix rightDirection bottomHalf topHalf
3dZcat -prefix wrong Direction topHalf bottomHalf

RightDirection

WrongDirection

A more useful application of 3dZcutup and 3dZcutup is during the stage of 3dDeconvolve, where each slice (or group of slices) can be run through 3dDeconvolve, and then stacked together to create the complete statistical dataset (the following is copied from the help file of 3dZcutup, since it is the better than any example I could come up with):

  foreach sl ( `count -dig 2 0 20` )
    3dZcutup -prefix zcut${sl} -keep $sl $sl epi07+orig

    # Analyze this slice with 3dDeconvolve separately

    3dDeconvolve -input zcut${sl}+orig.HEAD            \
                 -num_stimts 3                         \
                 -stim_file 1 ann_response_07.1D       \
                 -stim_file 2 antiann_response_07.1D   \
                 -stim_file 3 righthand_response_07.1D \
                 -stim_label 1 annulus                 \
                 -stim_label 2 antiann                 \
                 -stim_label 3 motor                   \
                 -stim_minlag 1 0  -stim_maxlag 1 0    \
                 -stim_minlag 2 0  -stim_maxlag 2 0    \
                 -stim_minlag 3 0  -stim_maxlag 3 0    \
                 -fitts zcut${sl}_fitts                \
                 -fout -bucket zcut${sl}_stats
  end

  # Assemble slicewise outputs into final datasets

  time 3dZcat -verb -prefix zc07a_fitts zcut??_fitts+orig.HEAD
  time 3dZcat -verb -prefix zc07a_stats zcut??_stats+orig.HEAD


What this will do is loop over twenty slices and perform 3dDeconvolve on each slice separately, and then reassemble both the fitts and stats datasets from all of the individual slices after they have been analyzed. This can help when the dataset is either extremely large, or your computer has relatively little memory.


Thanks to alert reader Landoska, who once cut his FMRI data into four slices instead of eight, because he wasn't hungry enough for eight slices. (rimshot)


Design Efficiency

An important consideration for any fMRI experimental design is how efficiently the BOLD response will be estimated, given the timing of events and jittering between trials of interest. Therefore, before doing any scanning, it is useful to do a dry run of the experimental task on a computer outside of the scanner, and then feed this timing information into the software analysis package that you will use. For example, in SPM this would involve Specifying the 1st-Level through the SPM GUI, but not estimating it.

In AFNI, this would involve using 3dDeconvolve with the -nodata command, telling AFNI that there is no functional data to read.

#!/bin/tcsh

3dDeconvolve -nodata 205 3                       \
    -num_stimts 18                                                         \
    -stim_times 1 stimuli/reg1.txt 'GAM'                               \
    -stim_label 1 reg1                                                 \
    -stim_times 2 stimuli/reg2.txt 'GAM'                          \
    -stim_label 2 reg2                                            \
    -stim_times 3 stimuli/reg3.txt 'GAM'                     \
    -stim_label 3 reg3                                       \
    -stim_times 4 stimuli/reg4.txt 'GAM'                          \
    -stim_label 4 reg4                                            \
    -stim_times 5 stimuli/reg5.txt 'GAM'                     \
    -stim_label 5 reg5                                       \
    -stim_times 6 stimuli/reg6.txt 'GAM'                           \
    -stim_label 6 reg6                                            \
    -stim_times 7 stimuli/reg7.txt 'GAM'                      \
    -stim_label 7 reg7                                       \
    -stim_times 8 stimuli/reg8.txt 'GAM'                          \
    -stim_label 8 reg8                                            \
    -stim_times 9 stimuli/reg9.txt 'GAM'                     \
    -stim_label 9 reg9                                       \
    -stim_times 10 stimuli/reg10.txt 'GAM'                            \
    -stim_label 10 reg10                                              \
    -stim_times 11 stimuli/reg11.txt 'GAM'                       \
    -stim_label 11 reg11                                         \
    -stim_times 12 stimuli/reg12.txt 'GAM'                                  \
    -stim_label 12 reg12                                                    \
    -stim_times 13 stimuli/reg13.txt 'GAM'                              \
    -stim_label 13 reg13                                               \
    -stim_times 14 stimuli/reg14.txt 'GAM'                       \
    -stim_label 14 reg14                                         \
    -stim_times 15 stimuli/reg15.txt 'GAM'                       \
    -stim_label 15 reg15                                        \
    -stim_times 16 stimuli/reg16.txt 'GAM'                       \
    -stim_label 16 reg16                                         \
    -stim_times 17 stimuli/reg17.txt 'GAM'                           \
    -stim_label 17 reg17                                             \
    -stim_times 18 stimuli/reg18.txt 'GAM'                    \
    -stim_label 18 reg18                                      \
    -x1D X.xmat.1D -xjpeg X.jpg                                \
    -x1D_uncensored X.nocensor.xmat.1D                                     \


(A couple notes about the above code: 3dDeconvolve -nodata should be followed by the number of time points, and the TR used. In this example, I had 205 TRs, with a TR of 3 seconds, for a total of 615 seconds per run. Also, the names of the regressors have been changed to generic labels.)

The output of this command will produce estimates of how efficiently 3dDeconvolve is able to fit a model to the data provided. The Signal+Baseline matrix condition is especially important, as any singularities in this matrix (e.g., perfect or very high correlations between regressors) will make any solution impossible, as this implies that there are potentially infinite answers to estimating the beta weights for each condition. You should see either GOOD or VERY GOOD for each matrix condition; anything less than that will require using the -GOFORIT command, which can override matrix warnings. However, this should be done with caution, and only if you know what you are doing. Which I don't, most of the time.

Also notice that a number labeled the "norm. std. dev." is calculated, which is the root mean square of the measurement error variance for each condition (not to be confused with plain variance). Higher unexplained variance is less desirable, and apparently this design results in very high unexplained variance. More details about how to mitigate this, as well as a more thorough discussion of measurement error variance in general, can be found in this AFNI tutorial.


 
This command also produces an X.xmat.1D file, which can then be used to either visualize the regressors via 1dplot, or to examine correlations between regressors. This is done with 1d_tool.py, which is the default in the uber_subject.py stream. 1d_tool.py will set the cutoff at 0.4 before it produces any warnings, although this can be changed with the -cormat_cutoff option.

In this example, I compared two designs, one using an exponential distribution of jitters (with more jitters at shorter durations), and another which sampled equally from jitters of different durations.While the default for 1d_tool.py is to show only correlations between regressors at 0.4 or above, this can be changed using the -cormat_cutoff command, as shown below:

1d_tool.py -cormat_cutoff 0.1 -show_cormat_warnings -infile X.xmat.1D

This command produced the following output for each design:

Equal Jitters

Exponential Jitters

Note that while both do not have correlations above 0.4, there are substantial differences when looking at the range of correlations between 0.1 and 0.3. In both cases, however, each design does a relatively good job of reducing correlations between the regressors. This should also be compared against the measurement error variance as described above, as a design may produce higher correlations, but lower measurement error variance.

Final note: There are programs to create efficient designs, such as Doug Greve's optseq and AFNI's RSFgen. However, this will produce designs that are expected to be used for each subject, which potentially could lead to ordering effects. Although the likelihood of this having any significant effect on your results is small, it is still possible. In any case, all of these things should be considered together when designing a scanning study, as it takes far less time to go through these steps than to end up with an inefficient design.

Duration Modulation in AFNI

Modulating signal by duration - usually reaction time (RT) - is increasingly common in fMRI data analysis, particularly in light of recent studies examining how partialing out RT can reduce or even eliminate effects in certain regions of the brain (e.g., anterior cingulate; see Grinband et al, 2010; Yarkoni et al, 2009). In light of these findings, it appears as though parametrically modulating regressors by RT, or the duration of a given condition, is an important factor in any analysis where this data is available.

When performing this type of analysis, therefore, it is important to know how your analysis software processes duration modulation data. With AFNI, the default behavior of their duration modulation basis function (dmBLOCK) used to scale everything to 1, no matter how long the trial lasted. This may be useful for comparison to other conditions which have also been scaled the same way, but is not an appropriate assumption for conditions lasting only a couple of seconds or less. The BOLD response tends to saturate over time when exposed to the same stimulation for an extended period (e.g., block designs repeatedly presenting visual or auditory stimulation), and so it is reasonable to assume that trials lasting only a few hundred milliseconds will have less of a ramping up effect in the BOLD response than trials lasting for several seconds.

The following simulations were generated with AFNI's 3dDeconvolve using a variation of the following command:

3dDeconvolve -nodata 350 1 -polort -1 -num_stimts 1 -stim_times_AM1 q.1D 'dmBLOCK' -x1D stdout: | 1dplot -stdin -thick -thick

Where the file "q.1D" contains the following:

10:1 40:2 70:3 100:4 130:5 160:6 190:7 220:8 250:9 280:30

In AFNI syntax, this means that event 1 started 10 seconds into the experiment, with a duration of 1 second; event 2 started 40 seconds into the experiment with a duration of 2 seconds, and so on.

The 3dDeconvolve command above is a good way to generate simulation data, through the "-nodata" option which tells 3dDeconvolve that there is no functional data to process. The command tells 3dDeconvolve to use dmBLOCK as a basis function, convolving each event with a boxcar function the length of the specified duration.

Running this command as is generates the following graph:

As is expected, trials that are shorter are scaled less, while trials lasting longer are scaled more, with a saturation effect occurring around 8-9 seconds.

Running 3dDeconvolve with a basis function scaling the signal change in each to 1 is done with the following:

3dDeconvolve -nodata 350 1 -polort -1 -num_stimts 1 -stim_times_AM1 q.1D 'dmBLOCK(1)' -x1D stdout: | 1dplot -stdin -thick -thick

And generates the following output:



Likewise, the ceiling on the basis function can be set to any arbitrary number, e.g.:

3dDeconvolve -nodata 350 1 -polort -1 -num_stimts 1 -stim_times_AM1 q.1D 'dmBLOCK(10)' -x1D stdout: | 1dplot -stdin -thick -thick


However, the default behavior of AFNI is to scale events differently based on different duration (and functions identically to the basis function dmBLOCK(0)). This type of "tophat" function makes sense, because unlimited signal increase as duration also increases would lead to more and more bloodflow to the brain, which, taken to its logical conclusion, would mean that if you showed someone flashing checkerboards for half an hour straight their head would explode.

As always, it is important to look at your data to see how well your model fits the timecourse of activity in certain areas. While it is reasonable to think that dmBLOCK(0) is the most appropriate basis function to use for duration-related trials, this may not always be the case.

These last two figures show the same subject analyzed with both dmBLOCK(0) and dmBLOCK(1). The underlying beta values for each do not differ significantly, although there is some variability in how much they differ in distinct cortical areas, and small but consistent changes in variability can lead to relatively large effects at the second level.

The image on the left hasn't been masked, but the underlying beta estimates should be the same in either case.

dmBLOCK(0)
dmBLOCK(1)