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 neuronal response (physiological effect) at 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 one specific regressor in the model, a second-order regressor of the interaction.

Suppose we are interested in finding out the connectivity of an ROI for a contrast between two conditions, A and B. MyInput+orig is the input file of a subject with which you have already run regression analysis with 3dDeconvolve. If there are multiple runs, run steps (1)-(4) for each run separately. The seed region can be an ROI, a peak voxel or a sphere around the peak voxel. The procedure is the same regardless of the experiment type, event-related or block design.

1. For each run extract the average time series of the ROI

3dmaskave -mask ROI -MyInput+orig > Seed.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 MyInput+orig > Seed.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'.

2. Remove the trend from the seed time series

3dDetrend -polort ? -prefix SeedR Seed.1D

Note: 3dDetrend only takes rows as input, so if you have an input file with a column in Seed.1D, add \' to "Seed.1D":

3dDetrend -polort ? -prefix SeedR Seed.1D\'

The output SeedR.1D is a one-row text file. Convert the one-row time series to one column:

1dtranspose SeedR.1D Seed_ts.1D

[Note - Instead of doing steps 1 and 2 for each run separately, an alternative is to reverse the order of the two steps: first use 3dSynthesize to extract the drift effect (baseline part), and subtract the trend from the original time series, and then extract the ROI from the concatenated time series. You may still need to break this ROI time series into multiple runs for later processing. Actually I prefer this approach because it accounts for other effects more properly in addition to the trend.]

3. Run deconvolution on the seed time series

First generate the impulse response function:

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

Then run:

3dTfitter -RHS Seed_ts.1D -FALTUNG GammaHR.1D Seed_Neur 012 0

Plot out Seed_Neur.1D and see if it looks reasonable comparing to the stimulus presentation in the experiment. You may want to experiment with different penality functions (check details of option -FALTUNG in 3dTffitter -help).

4. Obtain the interaction regressor

First create a 1D file, AvsBcoding.1D,  with 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). If you only consider one condition A, there are two options: (1) If the baseline condition more or less matches up with condition A, make it a constrast with condition A (coded with 1's) versus baseline (coded with -1's); (2) If you don't believe that there is modulation (interaction) between the seed and the baseline, code condition with 1's and all other time points with 0's.

1deval -a Seed_Neur.1D\' -b AvsBcoding.1D -expr 'a*b' > Inter_neu.1D

The interaction is created as

waver -GAM -peak 1 -TR ?  -input Inter_neu.1D -numout #TRs > Inter_ts.1D

5. Concate the regressors if there are multilple runs

Run cat separately on Seed_ts.1D (final output from step(2)) and Inter_ts.1D (final output from step(4)), and use the 2 concatenated 1D files for the next step.

6. Regression analysis

Basically add two more regressors from step (5) to your original 3dDeconvolve script by using option -stim_file, and add option -rout in 3dDeconvolve since the correlation coefficient for regressor Inter_ts.1D will be taken for group analysis. That is, include all of the original regressors (of interest and no interest) plus the two new ones. That way all sources of variabilities in the data would be properly accounted for.

7. Convert the correction coefficients for interaction to Z scores through Fisher transformation

3dDeconvolve can only output coefficient of determination R², not correlation coefficient R itself. So 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+orig'[SubbrickForR²]' -b Corr_subj01+orig'[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+orig -expr 'log((1+a)/(1-a))/2' -prefix Corr_subj01_Z

8. Repeat the above steps for all other subjects

9. Convert the brains of Z values to common space (e.g. tlrc)

10. Run group analysis with 3dttest on the Z values of the interaction effect: Keep in mind a high t value in the output of 3dtttest only indicates r is with some significance level different from 0, does not necessarily mean a high correlation.And if you want the correlation coefficient at group level, you can convert the z-score (sub-brick #0 in the output from the group analysis) back to r (the reverse conversion of step 9 above):