HRF estimation with 3dMSS

paulgeha · August 16, 2024, 5:52pm

Dear AFNI expert,

I have task based BOLD data where we delivered two different tastes (ms and tls) as well as a rinse (rns). The ms and tls had 4 blocks each per run with the following onsets and durations in seconds:
ms

0:50	
62.029: 36	
176.814: 73	
412.481: 54

and
tls

and finally the rinse rns has these onsets with 3 second durations: (49.9, 98.2,164.7,249.3, 315.8,400.4, 466.9, and 515.2 seconds). I also 17 different nuisance regressors including 6 motion regressors and their derivatives and 5 components from the CSF.

My data is already preprocessed in FSL and I will would like to use the preprocessed data to do an HRF estimation per the Chen et al., 2023 paper; I am comparing the response differences (ms minus tls) between patients and control. Are the following steps at the single subject level sensible before going to 3dMSS?

Below are the codes I used for one subject:
step1 :

3dDeconvolve -input denoised_func_data_nonaggr.nii.gz \
    -mask brain_mask.nii.gz \
    -polort A \
    -num_stimts 26 \
    -stim_times 1 msA_block1.1D 'TENT(0, 50, 51)' -stim_label 1 HRF_MsA_Block1 \
    -stim_times 2 msA_block2.1D 'TENT(0, 36, 37)' -stim_label 2 HRF_MsA_Block2 \
    -stim_times 3 msA_block3.1D 'TENT(0, 73, 74)' -stim_label 3 HRF_MsA_Block3 \
    -stim_times 4 msA_block4.1D 'TENT(0, 54, 55)' -stim_label 4 HRF_MsA_Block4 \
    -stim_times 5 tlsA_block1.1D 'TENT(0, 54, 55)' -stim_label 5 HRF_TlsA_Block1 \
    -stim_times 6 tlsA_block2.1D 'TENT(0, 54, 55)' -stim_label 6 HRF_TlsA_Block2 \
    -stim_times 7 tlsA_block3.1D 'TENT(0, 73, 74)' -stim_label 7 HRF_TlsA_Block3 \
    -stim_times 8 tlsA_block4.1D 'TENT(0, 36, 37)' -stim_label 8 HRF_TlsA_Block4 \
    -stim_times 9 rnsA1.1D 'BLOCK(3,1)' -stim_label 9 RNS \
    $(for i in $(seq 0 16); do echo "-stim_file $(($i + 10)) mc_csf_comcorr.txt'[$i]' -stim_base $(($i + 10)) -stim_label $(($i + 10)) Nuisance$(($i + 1)) "; done) \
    -fout -tout -x1D X.xmat.1D -xjpeg X.jpg \
    -x1D_stop

and then Step 2

3dREMLfit -matrix X.xmat.1D \
	-input denoised_func_data_nonaggr.nii.gz \
	-mask brain_mask.nii.gz \
	-Rbuck results_REML.nii \
	-Rvar results_REMLvar.nii \
	-verb \
	-GOFORIT 3

Thanks a lot for your help!

Paul Geha

Gang · August 17, 2024, 9:41pm

Paul,

I’d like to clarify a few points to ensure I fully understander the context:

Temporal Resolution: What is the temporal resolution (i.e., TR) of your fMRI data?
Hemodynamic Response Estimation: If feasible, estimating the hemodynamic response through a data-driven approach is generally more reliable than assuming a fixed, potentially inaccurate HRF. You could even use the tent basis function instead of relying on BLOCK(3,1) for the rinse phase. Given that the two task phases 'ms' and 'tls' have different durations, how do you plan to compare their estimated hemodynamic responses?
Covariate Handling: It seems you’ve preemptively removed 17 covariates related to head motion effects and CSF signals. I would recommend directly incorporating these covariates into the individual-level regression model rather than removing them beforehand. This is because pre-removal of these covariates could distort or bias your parameter estimates. In addition, two other issues are worth pondering:
- Task-Induced Motion: If head motion is partly induced by the task conditions, it can complicate and potentially bias the parameter estimation. Incorporating estimated motion effects as covariates might not be as effective as anticipated and could inadvertently distort the modeling process.
- CSF-Related Signals: Some evidence suggests that CSF-related signals may contain task-induced BOLD responses, so I’d be cautious about including these as covariates.

Gang Chen

paulgeha · August 17, 2024, 11:33pm

Hi Gang,

thanks for your reply.

My TR is 1 second.
I will add tent basis function for the rinse phase. However, I am actually not sure how to compare ms and tls; this is where I need advice I guess. As per your paper I suspect a difference between patients and controls in HRF.
I am new to AFNI so I thought that the removal of the covariates in the code I wrote is integrated into the regression that also includes ms, tls, and rns. It seems it is not. How would I put them all in one model. Thanks for that tip on the CSF. I will run the analysis without CSF signals ; I did see a correlation between the global signal (which I am not removing) and the CSF component.
Thanks a lot!

Paul

Gang · August 19, 2024, 1:05pm

Paul,

I am actually not sure how to compare ms and tls; this is where I need advice I guess. As per your paper I suspect a difference between patients and controls in HRF.

Do all participants have the same trial durations for each of the two conditions ('ms' and 'tls')? Since you're estimating the hemodynamic response for each trial individually, are you planning to compare patients and controls at the trial level or condition level?

-GOFORIT 3

Regarding the use of the '-GOFORIT 3' option, could you clarify the root cause of the collinearity concern you're addressing?

Gang Chen

paulgeha · August 19, 2024, 1:17pm

hi Gang,

Yes all subjects have the same trial durations and same number of blocks.

I used GOFORIT because my ms and tls are collinear unfortunately. Please see the uploaded design from FSL.
Thanks
Paul

Gang · August 19, 2024, 1:40pm

Paul,

If you don’t add -GOFORIT 3, does 3dDeconvolve or 3dREMLfit give you any warnings about collinearity?

The challenge lies in the differing trial-level durations of the 'ms' and 'tls' conditions. At the individual level, how do the estimated HDRs appear? One approach is to assess the magnitude of the estimated HDRs at the trial level and then compare the two conditions at the group level based on these trial-level magnitudes.

Gang Chen

paulgeha · August 19, 2024, 2:28pm

Gang,
yes, if I do not use GOFORIT I get a collinearity warning.

I am not sure I exactly follow what you mean: how do I estimate HDR in this case? do I use a canonical HRF and compare the groups?

Thanks again!
Paul

Gang · August 19, 2024, 3:13pm

I assume that rinsing occurs at the end of each trial. If that’s the case, it becomes difficult to reliably estimate the impact of rinsing because the BOLD responses from each trial and the rinsing phase overlap, leading to collinearity issues. Simply using the -GOFORIT 3 option won’t resolve this problem; it may just obscure it.

If you choose to estimate the trial-level BOLD response directly, I recommend against separately modeling the rinsing phase. Instead, append the rinsing phase to each trial and adjust your model specification accordingly—for instance, changing TENT(0, 50, 51) to TENT(0, 65, 66). Then, examine the morphology of the estimated HDRs to determine whether the magnitude of the HDR profile could be used for group-level analysis.

Alternatively, you could use the canonical approach as you mentioned, with the option dmUBLOCK in 3dDeconvolve (see here). Even then, you may still need to merge the rinsing phase with each trial.

You might try both modeling approaches and compare the consistency of the results at both the individual and group levels.

Gang Chen

paulgeha · August 19, 2024, 5:37pm

Gang,

I will model the rinse with the ms and tls but I think the collinearity is also coming from ms and tls correlation unfortunately.

Just so that I understand well; when you say "examine the morphology of the estimated HDR". Do you mean using 3dMSS?

Thanks
Paul

Gang · August 19, 2024, 7:25pm

I think the collinearity is also coming from ms and tls correlation unfortunately.

Unfortunately, this seems to be a fundamental experimental design issue. It's crucial to randomize and jitter the inter-trial intervals, ideally across conditions as well. A model cannot fix this issue retrospectively. If this is a pilot study, consider improving the design for the formal study. With the current data, one possible compromise is to shorten the recovery phase using something like TENT(0, 55, 56) and assess its effectiveness.

Just so that I understand well; when you say "examine the morphology of the estimated HDR". Do you mean using 3dMSS?

Given the variability in trial durations across trials and conditions, directly comparing HDRs based on their profiles using 3dMSS becomes challenging. Instead, consider extracting or summarizing the entire HDR profile by its magnitude to simplify group-level comparisons (e.g., using a simple model like two-sample t-test).

Gang Chen

paulgeha · August 19, 2024, 7:54pm

Ok thanks Gang!

paulgeha · August 19, 2024, 11:35pm

Gang,
one last note; my ms and tls tasks had 4 blocks each as I mentioned and are all the same for subjects who received the same protocol. We however have 4 different protocols A,B,C, and D and subjects received a random combination of A and B or C and D. The block durations or numbers do not differ betwen A-D but the inter-block separation is totally different.

I think this is what might partially answer your concerns.

Best

Paul

Gang · August 20, 2024, 1:40pm

Thanks for the clarification. Do the four blocks for each task condition correspond to four different protocols? Are you planning to compare the two groups for each protocol separately at the group level?

Another approach for group-level analysis is to focus on the major segment of the estimated HDRs. For example, you could extract the time points from, example, 0 to 36 seconds from the individual-level HDRs and compare that segment at the group level using 3dMSS.

Gang Chen

paulgeha · August 20, 2024, 1:51pm

Gang,

no we plan to combine the protocols; originally we planned to do the simple unpaired t-test (ms-tls)pat -(ms-tls)cont. The subject level contrast can come from A and B or C and D. One more piece of information. Since these were taste delivery blocks the delivery was not continuous within a block but it was through 2 pumps delivering the tastant over a 9 s period with < 1 s short inter-delivery interval within a block.

Once I figure out the lower level analysis with AFNI I can try the 3dMSS.

Thanks again!
Paul Geha.

Gang · August 20, 2024, 1:58pm

Paul,

Since these were taste delivery blocks the delivery was not continuous within a block but it was through 2 pumps delivering the tastant over a 9 s period with < 1 s short inter-delivery interval within a block.

In that case, it’s even more crucial to estimate the HDRs rather than blindly relying on a canonical HRF. By estimating the HDRs, you can visually inspect them and determine which segment you want to focus on at the group level.

Gang Chen

paulgeha · August 21, 2024, 5:21pm

Gang,
so to estimate the HDR, should I just try it on few subjects using 3dDeconvolve and then apply it to all of them ?

Thanks
Paul

Gang · August 21, 2024, 8:03pm

That sounds reasonable to me.

Gang Chen