How to use HRF with shorter latency for deconvolving in afni?

AFNI Message Board
1

Hi afni experts，
I am interested in the superior colliculus (sc) and sc has an optimised HRF with a 4s peak (Wall, M. B., Walker, R., & Smith, A. T. (2009). Functional imaging of the human superior colliculus: an optimised approach. Neuroimage, 47(4), 1620-1627.). The optimised HRF used in the reference is the standard HRF which is a synthetic function produced by combining two gamma functions (Friston et al., 1998; Handwerker
et al., 2004) and is implemented in BrainVoyager QX.
I don’t know how to use the same combined two gamma function with a 4s peak as HRF in afni. Though I notice the ‘GAM(p,q)’ option of Rmodel in 3dDeconvolve and the latency can be changed by changing the value of p*q. But the expression used in this option is just a single gamma function.
So :

  1. How to use the same synthetic function with a 4s peak as HRF in afni?
  2. If it is too complicated to use the synthetic function mentioned in the reference, ‘GAM(p,q)’ can also be accepted. But how to set the two parameters ‘p’ and ‘q’ properly? Could ‘p’ and ‘q’ be set arbitrarily to satisfy the equation p*q=4 if I need a 4s peak HRF? For example, p=2,q=2; p=8.6,q=0.4651 or p=7.3126,q=0.547. Are these combinations of ‘p’ and ‘q’ all OK? What do the two parameters (p & q) model in the gamma function? The magnitudes of the peak, time-to-onset or something else?
    Thank you in advance!
2

Unfortunately, the article referenced does not give the formula for the HRF they used, they just say it is implemented in BrainVoyager. So it is a little difficult to fully address your question.

There are a couple ways you can proceed, from easy to hard.

FIRST method: GAM(p,q) is described in detail in the output of 3dDeconvolve -help:
https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dDeconvolve.html


     'GAM(p,q)'    = 1 parameter gamma variate                         
                               (t/(p*q))^p * exp(p-t/q)                      
                             Defaults: p=8.6 q=0.547 if only 'GAM' is used   
                          ** The peak of 'GAM(p,q)' is at time p*q after    
                             the stimulus.  The FWHM is about 2.3*sqrt(p)*q.
                   ==> ** If you add a third argument 'd', then the GAM  
                              function is convolved with a square wave of    
                              duration 'd' seconds; for example:             
                                'GAM(8.6,.547,17)'                           
                              for a 17 second stimulus.

So if you want a peak to occur at D s, then pq=D should be chosen. If you want a (full) width of W s, then 2.3sqrt(p)q=W should be chosen. So D/W=sqrt(p)/2.3 or p=(2.3D/W)^2 and then q=D/p. For example, D=4 and W=6 gives p=2.35 q=1.70. You can visualize this function in AFNI commands via

1deval -num 200 -del 0.1 -expr 't^2.35*exp(-t/1.7)' | 1dplot -stdin -del 0.1

If you have a stimulus of duration 3 s (as in the cited paper), then using ‘GAM(2.35,1.7,3)’ would be appropriate for this model.

The drawback of this method, relative to the method of the cited paper, is that there is only one gamma variate (please do NOT call these gamma functions, as is commonly done, since those are something else entirely), and if you want to duplicate the analysis of the paper, you need two such functions added together appropriately.

At this point, I’m tired of typing, so will respond to this posting later when I gather my thoughts and rest my fingers.

3

The following changes are in the AFNI source code now, but will not be available in binaries until Tuesday morning 09 Jan 2018 (Tuesday evening in China).

I had planned to explain how to specify a two gamma variate model with the ‘EXPR’ model in 3dDeconvolve. But this is complicated and easy to use incorrectly. So instead, I modified the program to have a new model:

TWOGAM(p1,q1,r,p2,q2,d)

where p1,q1 are the parameters for the first GAM function; p2,q2 are the parameters for the second GAM function, r is the coefficient applied to the second one, and d is the duration (if d is not given, d is 0) as in

TWOGAM = GAM(p1,q1,d) - r * GAM(p2,q2,d)

Notice that parameter r is used to subtract the second GAM function (i.e., models undershoot).

Also, the GAM and TWOGAM functions have new variants that allow you to specify the time-to-peak (K) and the full-width-at-half maximum (W) directly, as in

TWOGAMpw(K1,W1,r,K2,W2,d)

Here, the ‘pw’ means ‘peak and width’. A reasonable looking model (with 0 duration) would be

TWOGAMpw(3,6,0.2,10,12)

You can use 3dDeconvolve to directly plot this, with a command like:


3dDeconvolve -num_stimts 1 -polort -1 -nodata 81 0.5               \
                     -stim_times 1 '1D: 0' 'TWOGAMpw(3,6,0.2,10,12)' \
                     -x1D stdout: | 1dplot -stdin -THICK -del 0.5

Using 3dDeconvolve and 1dplot together like this will let you adjust the parameters to get a response model that fits your needs and vision.

4

Hi Bob,
I am very very appreciated for what you did！

I don’t know how these parameters of the two gamma function specified in SPM or BrainVoyager, so I have some more questions:

  1. Is the variant d the time duration of stimuli or not? If each event lasts 2s in my event-related experiment, the d should be set to 2?
  2. To set these parameters of TWOGAMpw correctly, the signal change of my data is calculated firstly. I find the averaged signal change has a positive peak at the 4th second from stimuli onset and the width corresponding to this peak is 10s (from 0s to 10s from the stimuli onset). The averaged signal change also have a undershoot with a peak at the 12th second from stimuli onset and width of 4s (from 10s to 14s from stimuli onset) not including the part shared with the first peak.
    So the TWOGAMpw could be set like this:TWOGAMpw(4,5,r,12,7,2)?
  3. How to calculate the r of the undershoot from the signal change?

Thank you very much again!

5

If you want to duplicate the work in the paper by Wall (et al.), then you have to find out the parameters that BrainVoyager uses. I cannot find that information after spending all of 2 minutes with Google.

You can choose the K and W (peak and width) parameters based on looking at data, or at the models Wall used. Probably the exact values won’t matter much in the final results. The r value determines the size of the undershoot, and is usually in the range 0.1 to 0.2. Again, the exact value probably does not matter much.

You can plot your sample waveform with a command like this (here I’ve set r = 0.2):


3dDeconvolve -num_stimts 1 -polort -1 -nodata 51 0.5 \
   -stim_times 1 '1D: 0' 'TWOGAMpw(4,5,0.2,12,7)'       \
   -x1D stdout: | 1dplot -stdin -THICK -del 0.5

Making plots like this will let you decide quickly if the parameters make sense to you. The synthetic TR used is 0.5, to make the graphs look smooth.

6

Thank you, Bob!
I notice that the parameter d is omitted in your code (TWOGAMpw(4,5,0.2,12,7)). Do you think the duration should be set to 0 in my case though the stimuli duration is 2s?
By the way, the “TWOGAMpw” option is still unavailable although my afni binaries were updated using “@update.afni.binaries -defaults” this morning (at 18:00 on Wednesday evening 10 Jan 2018 in eastern American).

7

Thanks for letting us know. I guess you are using the
linux_openmp_64 package, which did not get built because
that machine was down for some reason. Those should be
ready in an hour or two.

  • rick
8

Please try updating again, the binaries should be ready.

  • rick