Typically subjects perform some tasks, or are under some conditions, or go through some stimuli during the scanning session. Thus different from simple correlation analysis, we may want to account for the effect of a specific task/condition/stimulus on the connectivity model in the analysis, and the interaction between this psychological effect (task/condition/stimulus) and the seed region. In other words, the integration process at mutliple levels in the brain is dependent on the context of tasks/conditions/stimuli. Other than the seed time course, such a context-dependent correlation analysis, aka phychophysiological interaction or physiophysiological interaction (PPI), requires the insertion of two specific regressors in the model, a first-order regressor of the hemodynamic response function (HRF) of the task/condition/stimulus and a second-order regressor of the interaction. Again an example below is only a suggestive process.

Sippose we are interested in finding out the connectivity of an ROI for a contrast between two conditions, A and B. If you only have one condition, make some adjustment.

(1) Suppose you have all the beta weights (option -cbucket in 3dDeconvolve) and design matrix (-x1D_uncensored in 3dDeconvolve) available. Extract the effect of each condition for every subject:

3dSynthesize -cbucket ... -matrix ... -select baseline RegressorsOtherThanAandB ... -prefix AnythingExceptCondAandCondB

3dcalc -a InputFor3Deconvolve -b AnythingExceptCondAandCondB -expr 'a-b' -prefix EffectsOfCondAandCondB

(2) Convert the final output from (1) to common space

adwarp -apar anat_subj01+tlrc -dpar 'EffectOfCondAandCondB_subj01+orig' -dxyz 3 -prefix EffectsOfCondAandCondB_subj01

(3) Extract the average time series of the ROI

3dmaskave -mask ROI -EffectOfCondAandCondB_subj01+tlrc > EffectOfCondAandCondB_subj01.1D

or get the time course of the voxel whose t value peaks within the ROI [at (46, -29, 21) in this case]

3dmaskdump -noijk -dbox 46 -29 21 EffectOfCondAandCondB_subj01+tlrc > EffectOfCondAandCondB_subj01.1D

Note:

A. -noijk avoids writing out 3 coordinate index numbers at the beginning to the output file.
B. Regarding -dbox, be careful with different coordinate system. Here is some detail about the issue. Also check out '3dmaskdump -help'.

Convert the one-row time series to one column:

1dtranspose EffectOfCondA_subj01_DT.1D SeedCondAandCondB_subj01.1D

(4) Regressors of no interest (optional)

Some people recommend that physiological data such as cardiac and respiratory rates be recorded during the scanning and be included as covariates in correlation analysis. If such physiological data are not available, the global signal might be a remedy:

3dAutomask -dilate 1 -prefix EPI_subj01_mask EPI_subj01+tlrc
3dmaskave -mask EPI_subj01_mask+tlrc EPI_subj01+tlrc > subj01_global.1D

(5) Obtain the contrast regressor

waver -GAM -peak 1 -TR 2  -input AvsBcoding.1D -numout #TRs > AvsB.1D

Here AvsBcoding.1D is a file of 0's (at those TR's where neither condition A nor B occurred), 1's (at those TR's where condition A occurred), and -1's(at those TR's where condition B occurred).

(6) Construction of interaction regressor

For a block design, the interaction between the seed region and the contrast is obtained by simply calculating the pair-wise product between the two time series (e.g., Heekeren et al., 2004):

1deval -a SeedCondA_subj01.1D -b AvsBcoding.1D -expr 'a*b' > InteractionAvsB.1D

Here AvsBcoding.1D is a file of 0's (at those TR's where neither condition A nor B occurred), 1's (at those TR's where condition A occurred), and -1's (at those TR's where condition B occurred). 

For an event-related design (Gitelman et al., 2003), we need to take the consideration of interaction at the neuronal level. Instead of obtaining the interaction regressor by simple multiplication at the BOLD signal level, we estimate the neuronal strength for the event through deconvolution, an opposite operation of convolution.

First generate the impulse response function:

waver -dt TR -GAM -inline 1@1 > GammaHR.1D

Then run deconvolution to get the response at neuronal level (you may want to change some parameters in the following command):

 3dTfitter -RHS SeedCondAandCondB_subj01.1D -FALTUNG GammaHR.1D CondAndCondBNeuronal 012 0.0

The interaction at neuronal level can be represented by the pair-wise product between the event timing file in 0s and 1s and the seed response, Neuron.1D,

1deval -a AvsBcoding.1D -b NeuronCondAndCondB -expr 'a*b' > NR_AvsB.1D

Again AvsBcoding.1D is a file of 0's (at those TR's where neither condition A nor B occurred), 1's (at those TR's where condition A occurred), and -1's(at those TR's where condition B occurred).

The effect due to this interaction at hemodynamic level is computed as

waver -GAM -peak 1 -TR 2  -input NR_AvsB.1D -numout #TRs > InteractionAvsB.1D

(7) Correlation analysis

3dDeconvolve -input EffectsOfCondAandCondB_subj01 \
-polort 0 \                                    # since drifting effect has already been removed in the input file ffectOfCondAandCondB_subj01
-num_stimts 4 \    # Maybe more than 4 in your case!!!
-stim_file 1 subj01_global.1D -stim_label 1 "RegressorOfNoInterest" \  # This is totally optional: global signal, physiological regressors, other regressors not considered for the connectivity analysis, etc.
-stim_file 2 EffectOfCondAandCondB_subj01.1D -stim_label 2 "CorrelatonWithSeed" \         # Correlation with the seed
-stim_file 3 AvsB.1D -stim_label 3 "AvsB" \                 # AvsB
-stim_file 4 InteractionAvsB.1D -stim_label 4 "Interaction" \              # Interaction between seed region and the contrast
-tout -rout -fitts fit_subj01 \
-bucket Corr_subj01

(8) Repeat the above steps for all other subjects

(9) Convert the correction coefficients for "Interaction" (regressor 6 above) to Z scores through Fisher transformation

3dDeconvolve can only output R², not correlation coefficient R itself. So if you want R, we need to take square root of R² and find out its sign based on the sign of its corresponding beta value:

3dcalc -a Corr_subj01+tlrc'[SubbrickForR²]' -b Corr_subj01+tlrc'[SubbrickForBeta]'-expr 'ispositive(b)*sqrt(a)-isnegative(b)*sqrt(a)'-prefix Corr_subj01R

Since correlation coefficients range from -1 to 1. To be able to run group analysis, Fisher's Z transformation formula can be used to reduce skewness and make the sampling distribution more normal when sample size is big enough: z =  (1/2) * ln((1+r)/(1-r)), where z is approximately normally distributed with mean r, and standard error 1/(n-3)^0.5 (n: sample size).

3dcalc -a Corr_subj01R+tlrc -expr 'log((1+a)/(1-a))/2' -prefix Corr_subj01_Z

(10) Group analysis on Z scores with 3dttet

References

Heekeren et al., A general mechanism for perceptual decision-making in the human brain. Nature Vol. 431, 14 October 2004.
Gitelman et al., Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution. NeuroImage, 19:200-207, 2003

* Back to Gang Chen's home page