Estimating the hemodynamic response (“Deconvolution”)

 

  1. Overview of the approach
  2. Tent basis functions
  3. Sine basis functions
  4. Contrasts and basis functions in AFNI
  5. Interactively exploring basis functions prior to analysis

1. Overview of the approach

[for a more mathematical discussion of this technique, see the 3dDeconvolve manual].

 

In the simple example discussed here, we used a Gamma variate function as a model of the shape of the hemodynamic response (HR), and simply used GLM to estimate the height, or amplitude of the response for different experimental conditions. This technique works well when the shape of the actual HR is close to that of a Gamma variate, which might be a reasonable assumption for simple, short stimuli, particularly in primary visual/auditory/motor cortex.

 

But in some experimental situations, we might expect the HR to have a markedly different shape, and even to vary across conditions, brain regions and/or subjects. In these situations, it is preferable to use a technique that estimates not only the amplitude of the HR, but also its shape. The process by which we try to infer the shape (and amplitude) of the HR for each experimental condition from the total MR signal can be considered a type of deconvolution. The ways in which deconvolution can be formed are numerous, but in AFNI (and in other fMRI analysis software) the method is to model the HR for each condition as a weighted linear sum of a set of basis functions.

 

To get an idea of how this works, it is worth looking at a few simplified examples. Say we have an actual HR that looks like the blue curve below:

 

 

We could model this HR with a weighted linear sum of the following 6 basis functions, B1 to B6 (each color represents one basis function):

 

 

We use GLM to fnd a set of weights, w1 to w6 such that:

 

w1B1 + w2B2 + w3B3 + w4B4 + w5B5 + w6B6 ~ HR

 

In this case, the best set of weights would be: -0.28, -0.72, 10.03, 8.91, 2.05, -1.48, which when used in the formula above give the estimated HR, shown with a broken yellow line in the top graph.

 

So the question poses itself: what type, and how many basis functions should we use to model a given expected HR? Remember that the purpose of using deconvolution in the first place is to allow flexibility to model different HR shapes. Accordingly, there are no simple answers to the question. Fortunately, AFNI provides a number of “ready made” basis function sets that do a reasonable job for a variety of possible HR shapes. Of these, we concentrate here on tent functions, and sine functions.

 

<back to top>

 

2. Tent Basis Functions.

The example that we just showed above used tent functions, which are simply “tent-shaped” functions that span the rnage of times over which we expect a given response to occur. AFNI allows us to specify the time span and number of tent functions used to model a given HR using the following syntax:

 

TENT(b,c,n)

 

Where b represents the start of the modeled response (relative to the stimulus onset), c represents the end of the modeled response (also relative to stimulus onset), and n gives the desired number of tent functions. Using our previous example, we might choose to model the reponse to positive images as:

 

-stim_times 1 pos_onsets.1D ‘TENT(0,20,8)’

 

This would thus allow us to model the response as occuring over the period from 0 to 20 seconds after each positive image onset, using 8 tent functions. An example is shown below, for a hypothesised HR in blue, and the fitted estimated HR in yellow. The 8 basis functions that were used to make the estimate are shown in the bottom panel.

 

 

<back to top>

 

3. Sine Basis Functions.

The other type of basis function that we will discuss here is the sine basis function. It turns out that it’s possible to model a variety of different HR shapes with a set of sine wave basis functions, each with a different frequency. This is done in AFNI using the following syntax:

 

SIN(b,c,n)

 

Where b represents the start of the modeled response (relative to the stimulus onset), c represents the end of the modeled response (also relative to stimulus onset), and n gives the desired number of tent functions. Using our previous example, we might choose to model the reponse to positive images as:

 

-stim_times 1 pos_onsets.1D ‘SIN(0,20,5)’

 

This would thus allow us to model the response as occuring over the period from 0 to 20 seconds after each positive image onset, using 5 sine functions. An example is shown below, for a hypothesised HR in blue, and the fitted estimated HR in yellow. The 5 basis functions that were used to make the estimate are shown in the bottom panel.

 

 

 

Unlike the tent function approach, each sine function spans the entire window of time over which the HR is modeled. This has some advantages; it is somewhat less susceptible to the influence of outliers in the data, which is particularly relevant for designs with a small number of trials per condition. Sine functions are also inherently smooth, and this imposes a constraint on the shpae of the HR which makes sense theoretically, since hemodynamic changes are also smooth. Finally, often fewer sine wave basis functions than tent basis functions are required to achieve a comparable degree of precision in the estimated HR. Having fewer basis functions means there are fewer parameters (weights) to estimate in the GLM, which improves the reliability of the estimates.

 

One caveat about using sine basis functions is that they can only model responses which start with 0 amplitude (i.e. the response has 0 amplitude at the first modeled time point). This is a reasonable assumption for many designs, but might not be for any particular design.

 

<back to top>

 

4. Contrasts and basis functions in AFNI

Normally when using a preconvolved (e.g. Gamma variate) model in individual subject GLM, we estimate a parameter that corresponds to the amplitude of the HR for each experimental condition. It is thus trivial to also estimate contrasts, for example the subtraction of one parameter estimate from another to give HR amplitude differences between conditions.

 

With basis function approaches to individual subject GLM, things are not so simple. When we estimate the HR using basis functions, the individual parameter estimates correspond to weighting factors, and have no clear intrinsic physical meaning. They are only meaningful when multiplied by their respective basis functions and added together. There is thus no meaningful way to calculate contrasts directly with the 3dDeconvolve command.

 

Fortunately, 3dDeconvolve does provide an option to output the estimated HR (i.e. the estimated weights multiplied by the basis functions and added together) for each condition and each voxel as a 3d+time AFNI brain image. Using the example above, where we specified the modeling of responses to positive images with:

 

-stim_times 1 pos_onsets.1D ‘SIN(0,20,5)’

 

we would also add the following line:

 

-iresp 1 pos_response

 

This would produce an AFNI dataset called pos_response+orig.HEAD/BRIK which would contain, for each voxel, a 20 second estimate of the HR to positive images. The time resolution of this dataset is by deault the same as the original data. To provide a 1-second time resolution, would could additionally add the following line:

 

-TR_times 1

 

The estimated HR datasets can then be viewed using the “graph” feature in AFNI.

 

Now that we have the estimated HR for each condition, we are in a position to calculate contrasts, which we can do using the AFNI program 3dcalc. But the question is, what do we want to contrast? If we wanted to compare the amplitudes of two responses, we might contrast the maxima of the two responses:

 

3dcalc –prefix pos_neg_ampl –a pos_response+orig –b neg_response+orig –expr “max(a)-max(b)”

 

Of course, this might not give very reliable estimates if the estimated HRs themselves aren’t very reliable, since we’re basing the calculation on single time points in each response. Another option is to decide apriori what range of HR time points are of interest (e.g. those likely to correspond to the peak of activation, based upon prior work or theroretical considerations), and then contrast the mean of those time points across conditions:

 

3dcalc –prefix pos_neg_TR2_TR4 –a pos_response+orig[2] –b pos_response+orig[3] –c pos_response+orig[4] \

-d neg_response+orig[2] -e neg_response+orig[3] -f neg_response+orig[4] –expr “mean(a,b,c)-mean(d,e,f)”

 

NOTE: Make sure you’re taking into account the time resolution (i.e. 1 second or 1 TR) of the HR datasets.

 

The type of contrast you decide to calculate is really dependent on your research question. The resulting contrasts can then be converted to percent signal change and registered to a standardised template for group analysis.

 

<back to top>

 

5. Interactively exploring basis functions prior to analysis

Regardless of which type of basis functions are chosen to model the HR for a given condition, a decision still needs to be made as to how many basis functions are required. As stated above, generally fewer sine basis functions will be required than tent basis functions. But the exact number to use depends on a few things. For one, the longer the expected HR, the more basis functions that will be required. Also, the more complex the expected HR, the more basis functions that are needed to model it.

 

It is always possible to specify lots of basis functions, and therefore in principle fit any arbitrary HR with high precision. However, this would not be a good idea. Remember that every basis function requires the estimate of an associated weight, and that the more weights that need to be specified, the less reliable any particular weight estimate will be. This will particularly be the case in fMRI experiments where the HRs are likely to overlap (e.g. fast event-related experiments), since in essence multiple weights will be estimated on the basis of the same set of data points. Also, by specifying lots of basis functions, it is more likely that the GLM will be able to come up with a good fit even for “noise” – the more that you can constrain your analysis on the basis of what you know are reasonable HR shapes, the better.

 

In summary, it is probably best to limit the number of basis functions to the minimum number required to reasonably accurately model the range of expected HR shapes in your experiment. To assist in this process, you can use the program BasisFunctions, which can be started from any Linux machine on the Keck network. Just type:

 

BasisFunctions

 

at the command prompt (please email Tom Johnstone if it does not appear to be working). This program allows you to specify hypothesised HRs, and how many tent or sine basis functions to estimate the HR with. The program then performs an estimation, plotting both the hypothesised HR and the fitted estimate. By adjusting the number and type of basis functions, and the shape of the HR, you should get a feel for the appropriate settings to use in 3dDeconvolve (HINT: you might also want to create some “noise” HRs, and see what happens to the estimated HR with different settings).

 

<back to top>