How to use 3dMSS to analyze task-fMRI data?

AFNI version info (afni -ver):
Precompiled binary macos_13_ARM_clang: Dec 7 2023 (Version AFNI_23.3.10 'Septimius Severus')

Hi experts and @ptaylor,

I recently tried to analyze my task-fMRI data with 3dMSS to extract more results in white matter. I have pre-processed the EPI and anatomy data in subject level by "afni_proc.py" and generated the beta effect image for each subject. I also performed group comparisons between the normal group and patients group by 3dttest++.

I have read the 3dMSS article (BOLD Response is more than just magnitude: Improving detection sensitivity through capturing hemodynamic profiles), example code (AFNI GitHub) and 3dMSS help(AFNI homepage). However, I am still not sure how to deal with my data with 3dMSS. Should I use 3dMSS process the data which has been pre-processed by afni_proc.py (including tshift align tlrc volreg mask blur scale regress blocks)? But according to my understanding, modeling the HRF with smooth splines should use the 4D data, which has been removed in the result of afni_proc.py. Or Should I only pre-process the data with the "tshift align tlrc volreg mask blur scale"blocks ? I am confused with the analysis flow with 3dMSS.

By the way, I found an 3dMSS option "-bounds", which was not used in "afni_proc.py". What is the function of this option? Should I use the same parameter [-2 2] as the example.

Thank you!

Best regard,
Roger

Hi Roger,

Should I use 3dMSS process the data which has been pre-processed by afni_proc.py (including tshift align tlrc volreg mask blur scale regress blocks)?

At the individual level, you can estimate the hemodynamic response at the condition level using basis functions such as TENT or CSPLIN in 3dDeconvolve/3dREMLfit. Then, take those regression coefficients (i.e., \beta's associated with those basis functions) as input for 3dMSS at the population level.

By the way, I found an 3dMSS option "-bounds", which was not used in "afni_proc.py". What is the function of this option? Should I use the same parameter [-2 2] as the example.

As explained in the 3dMSS help, the option -bounds is intended to remove those input values that are considered to be outliers. The specific bound values should be chosen based on the context.

Gang Chen

Hi Prof.Chen @Gang ,

Thank you for your detailed reply. I read the article "BOLD Response is more than just magnitude: Improving detection sensitivity through capturing hemodynamic profiles", and it is reported that smooth fitting (thin plate splines) had some advantages in detection sensitivity, reproducibility and compared to sampled HRF or TENT. So I supposed that it may need the 4D data.

However, it seems that it fits the data in population level, which I am curious about but can't understand. Could you provide a more detailed analysis protocol for task-fMRI data? Meanwhile, what is the difference between 3dDeconvolve and " afni_proc.py -regree block"? It seems that both of them could fit the 4D data with ideal HRF. How could I combine the preprocess (tshift align tlrc volreg blur mask scale) with 3dDeconvolve and 3dMSS? Thanks for your help.

Best,

Roger

Roger,

Could you provide a more detailed analysis protocol for task-fMRI data? Meanwhile, what is the difference between 3dDeconvolve and " afni_proc.py -regree block"? It seems that both of them could fit the 4D data with ideal HRF. How could I combine the preprocess (tshift align tlrc volreg blur mask scale) with 3dDeconvolve and 3dMSS?

To clarify, 3dMSS functions as a program that infers the population-level hemodynamic response. The following emphasis aims to clarifies the confusions:

  • The methodology in 3dMSS aligns with conventional one- and two-sample t-tests in fMRI experiments, with a crucial distinction. Unlike t-tests, which assume the hemodynamic response is represented by its magnitude (or a scalar relative to the canonical HRF) through the regression coefficient (\beta) at the individual level, 3dMSS requires the expression of the hemodynamic response through its overall shape, which is derived from each individual. This shape can be estimated using basis functions such as TENT or CSPLIN through 3dDeconvolve at the individual level.

  • To estimate the hemodynamic response individually, you can either directly specify the basis function using 3dDeconvolve with options like -stim_times ... TENT..., or convey the information through the processing script afni_proc.py ... -regress_basis TENT....

These points are also illustrated in this blog post. Please let us know if this explanation brings clarity.

Gang Chen

Hi Prof.Chen @Gang ,

Thank you for the information provided above. I have read the post and AFNI help about 3dDeconvolve. I analyzed my data with TENT in the individual level. Some of the code as follows:

afni_proc.py
....
        -regress_stim_times \
        /stimuli/ConditionA.txt   \
        //stimuli/ConditionB.txt     \
        -regress_stim_labels ConA ConB \
        -regress_basis 'TENT(0,10,5)'                    \
        -regress_opts_3dD  -jobs 8                                      \
                           -gltsym 'SYM: ConA - ConB'           \
                           -glt_label 1 A-B   
...

Then I checked the result and found each subject has 26 volumes, including one Full_Fstat, 11 volumes for ConA (including five Coef, five Tstat and Fstat), 11 volumes for ConB and three for A-B (Coef, Tstat and Fstat). But I am not sure what is my next step? I am going to try the following code with 3dMSS and dataTable.txt. If I want to compare the activation of ConA between the two group, could you help me check whether it is right? Meanwhile, if I want to compare the difference between ConA and ConB among the two groups, how should I write the code?

3dMSS -prefix $outputDir/MSS -jobs 10 \
      -mrr 's(age)+s(TR,by=grp)' \
      -qVars 'age, TR' \
      -mask $inputDir/mask_group.nii \
      -bounds -2 2 \
      -prediction $outputDir/@pred.txt \
      -dataTable $inputDir/@data.txt

data.txt:

subject age grp TR InputFile
s1 23 -1 0 stats.PA4.BRIK'[1]'
s1 23 -1 2.4 stats.PA4.BRIK'[2]'
s1 23 -1 4.8 stats.PA4.BRIK'[3]'
s1 23 -1 7.2 stats.PA4.BRIK'[4]'
s1 23 -1 9.6 stats.PA4.BRIK'[5]'
s2 24 -1 0 stats.PA5.BRIK'[1]'
s2 24 -1 2.4 stats.PA5.BRIK'[2]'
s2 24 -1 4.8 stats.PA5.BRIK'[3]'
s2 24 -1 7.2 stats.PA5.BRIK'[4]'
s2 24 -1 9.6 stats.PA5.BRIK'[5]'
....
s7 29 1 0 stats.PA10.BRIK'[1]'
s7 29 1 2.4 stats.PA10.BRIK'[2]'
s7 29 1 4.8 stats.PA10.BRIK'[3]'
s7 29 1 7.2 stats.PA10.BRIK'[4]'
s7 29 1 9.6 stats.PA10.BRIK'[5]'
....

pred.txt:

label   age    group   TR
s1 23 -1 0
s1 23 -1 0.3
s1 23 -1 0.6
s1 23 -1 0.9
s1 23 -1 1.2
s1 23 -1 1.5
s1 23 -1 1.8
...
s1 23 -1 9.6
s2 24 -1 0
s2 24 -1 0.3
...
s2 24 -1 9.6
...
...
s7 29 1 0
...
s7 29 1 9.6

However, when I executed these codes, I still got the error as "ERROR: The last column must be "InputFile" [without quotes] (or "Ausgang_val" in some cases)!!!". I tried changing the file names (which don't exist) in InputFile column, I got the same error again. So I don't know where is the problem.

By the way, is there any method to generate this table more convenient when I have lots of patients and time point?

Thank you for your help.

Best regards,

Roger

Hi Roger,

I have a question that whether this HRF profile analysis could be applied to block design tasks with an interval between each block.

While the BOLD response in a block design experiment is stronger, subtle variations in the response shape may still occur. For instance, the duration of the uprising phase may vary across regions, and the plateau phase may not remain flat. Therefore, HRF estimation can still enhance detection sensitivity.

to what degree of the randomization is needed for HRF profile analysis to reduce collinearity?

Regarding the degree of randomization needed for HRF profile analysis to reduce collinearity, randomization can help in two ways: by avoiding participants' anticipation and reducing collinearity among explanatory variables. It may suffice if these issues are mild in your case.

-regress_basis 'TENT(0,10,5)'

Assuming a TR of 2 s in your case, you might want to extend the HDR coverage. For example, if the block duration is 8 s, consider using 'TENT(0,22,12)'.

With two conditions A and B, you should have 13 estimated values for the contrast of A - B from the 3dDeconvolve output. Use these 13 sub-bricks as input for 3dMSS and fill them into your file data.txt .

is there any method to generate this table more convenient when I have lots of patients and time point?

As for generating tables more conveniently, there is currently no dedicated tool available. However, you can write a script in your preferred language to automate the generation of these tables.

Feel free to reach out if you have more questions.

Hope this helps,
Gang Chen