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.