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 p*q=D should be chosen. If you want a (full) width of W s, then 2.3*sqrt(p)*q=W should be chosen. So D/W=sqrt(p)/2.3 or p=(2.3*D/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.