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 multiple 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 psychophysiological 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.

The following illustrates the conventional analysis between the contrast of two conditions. However, a generalized approach may be better, and is discussed in the second half of this page.

Conventional approach

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.]

2a. If your stimulus onset times were not synchronized with TR grids, pick up a sub_TR, e.g., 0.1 seconds, replace the above 1dtranspose step and upsample seed time series by xx (original TR divided by sub_TR) times:

1dUpsample xx SeedR.1D\' > Seed_ts.1D

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).

3a. In case your stimulus onset times were not synchronized with TR grids, replace the above waver command with

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

4. Obtain the interaction regressor

First create a 1D (one column) 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, code condition with 1's and the baseline time point with -1'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

4a. If your stimulus onset times were not synchronized with TR grids, for each condition you can obtain the 1s (condition present) and 0s (condition absent) with the following command:

timing_tool.py -timing condition_timing_in_original_TR -tr sub_TR -stim_dur ... -run_len ... -min_frac ... -timing_to_1D ... -per_run -show_timing

And in the end you may need to down-sample the interaction time series back to TR grids by running

1dcat Inter_ts.1D'{0..$(xx)}' > Inter_ts_F.1D

where xx = original TR divided by sub_TR.

5. Concatenate the regressors if there are multiple 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. Group analysis

My limited experience showed that running group analysis with individual subject beta and the corresponding t-test with 3dttest/3dttest++/3dMEMA is more favorable than the following alternative approach with correlation coefficients. So in case you want to go with correlation coefficient, here are the steps:

a. 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

or more concisely,

3dcalc -a Corr_subj01R+orig -expr 'atanh(a)' -prefix Corr_subj01_Z

b. Repeat the above steps for all other subjects

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

d. 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):

3dcalc -a Grp_Result+tlrc'[0]' -expr '(exp(2*a)-1)/(exp(2*a)+1)' -prefix GrpR

or,

3dcalc -a Grp_Result+tlrc'[0]' -expr 'atanh(a)' -prefix GrpR

8. Interpretation of the effect of seed on a target region.

First of all, the PPI effect is independent of the typical effect (e.g., the contrast between conditions A and B). In other words, positive PPI effect has nothing to do with the sign of the contrast between conditions A and B. Secondly, a positive context-dependent (or interaction) effect means:

(1) the seed region tends to increase the contrast between the effects of A and B on the target region; or
(2) the contrast between the effects of A and B  tends to increase the effect of the seed region on the target region; or
(3) both (1) and (2).

Vice versa for the negative value.

A generalized approach

The generalized analysis method is discussed in the following two papers:

McLaren, D. G., Ries, M. L., Xu, G., & Johnson, S. C. (2012). A Generalized Form of Context-Dependent Psychophysiological Interactions (gPPI): A Comparison to Standard Approaches. NeuroImage, 1–41. doi:10.1016/j.neuroimage.2012.03.068

Cisler, J. M., Bush, K., & Steele, J. S. (2013). NeuroImage, 1–11. doi:10.1016/j.neuroimage.2013.09.018

Suppose there are multiple conditions (e.g. positive, negative, and neutral) or a factorial design (e.g., 9 conditions: 3 categories (human, animal, object) by 3 valences (positive, negative, and neutral)). We need to create an interaction regressor for each condition separately. Any contrasts among the conditions are not performed at the individual subject level, but at the group level. More specifically, steps 1 and 3 above are the same for all conditions. For each condition, follow steps 2, 4, and 5 above except that at step 4 there are no more -1's (since contrast is not considered at the individual level).

In the end, add all the interaction regressors (as many as the number of conditions of interest) plus the seed time series in the original individual subject regression script (which should include the regressors of average effect for the conditions in the model).

At the group level, you may run paired t-test, ANOVA, or even ANCOVA with the effect estimates associated with those interaction regressors, allowing for teasing apart the effects in more dimensions.