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.

From now on, we recommend a generalized analysis method discussed in the following two papers (the traditional approach is currently considered outdated):

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.

1. For each run extract the average time series of the ROI (if afni_proc.py was used before, consider using the pb04.* file from the pipeline as input here):

3dmaskave -mask ROI -quiet 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, conditionA.1D,  with 0's (at those TR's where condition A does not occur), 1's (at those TR's where condition A occurs). Unlike the traditional PPI, no contrasting is performed here (thus no -1's).

1deval -a Seed_Neur.1D\' -b conditionA.1D -expr 'a*b' > Inter_neuA.1D

The interaction is created as

waver -GAM -peak 1 -TR ?  -input Inter_neuA.1D -numout #TRs > Inter_A.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_A.1D'{0..$(xx)}' > Inter_Ad.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

Take beta value associated with each interaction regressor to group analysis (paired t-test, ANOVA, or even ANCOVA).

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.

Script template

The following is an example that demonstrates the pipeline for generalized PPI with three conditions at the individual level. Modify it to suite your own situation.

set subj = Sub1
# there are 5 runs
set nruns = 5
# number of time points per run in TR
set n_tp = 480
set TR = 2
# up-sample the data because of stimulus duration of 3s
set sub_TR = 1
# seed label
set sd = amygdala
# three conditions
set condList = (A B C)

# create Gamma impulse response function
waver -dt ${sub_TR} -GAM -peak 1 -inline 1@1 > GammaHR.1D

# for each run, extract seed time series, run deconvolution, and create interaction regressor
foreach cc (`count -digits 1 1 $nruns`)
   3dmaskave -mask ROI+orig -quiet pb04.$subj.r0${cc}.scale+tlrc > Seed${cc}${sd}.1D
   3dDetrend -polort 2 -prefix SeedR${cc}${sd} Seed${cc}${sd}.1D
   1dtranspose SeedR${cc}${sd}.1D Seed_ts${cc}${sd}D.1D
   rm -f SeedR${cc}${sd}.1D
   1dUpsample 2 Seed_ts${cc}${sd}D.1D > Seed_ts${cc}${sd}.1D
   3dTfitter -RHS Seed_ts${cc}${sd}.1D -FALTUNG GammaHR.1D
   Seed_Neur${cc}${sd} 012 -1
   foreach cond ($condList)
      head -${cc} stimuli/Allruns_stim_${cond}_time.1D |tail -1 > tmp.1D
      waver -dt ${sub_TR} -FILE ${sub_TR} one1.1D -tstim `cat tmp.1D` -numout ${n_tp} > ${cond}${cc}${sd}.1D
      rm -f tmp.1D
      1deval -a Seed_Neur${cc}${sd}.1D\' -b ${cond}${cc}${sd}.1D -expr 'a*b' > Inter_neu${cond}${cc}${sd}.1D
      waver -GAM -peak 1 -${TR} ${sub_TR} -input Inter_neu${cond}${cc}${sd}.1D -numout ${n_tp} > Inter_hrf${cond}${cc}${sd}.1D
      1dcat Inter_hrf${cond}${cc}${sd}.1D'{0..$(2)}' > Inter_ts${cond}${cc}${sd}.1D
   end
end

# catenate the two regressors across runs
cat Seed_ts?${sd}D.1D > Seed_ts${sd}.1D
cat Inter_ts${cond}?${sd}.1D > Inter_ts${cond}${sd}.1D

# re-run regression analysis by adding the two new regressors

3dDeconvolve -input pb04.$subj.r??.scale+tlrc.HEAD                      \
    -polort A                                                           \
    -mask full_mask.$subj+orig                                          \
    -num_stimts 13                                                      \
    -stim_times 1 stimuli/Allruns_stim_A_time.1D 'BLOCK(3,1)'     \
    -stim_label 1 A                                               \
    -stim_times 2 stimuli/Allruns_stim_B_time.1D 'BLOCK(3,1)'     \
    -stim_label 2 B                                               \
    -stim_times 3 stimuli/Allruns_stim_C_time.1D 'BLOCK(3,1)'     \
    -stim_label 3 C                                               \
    -stim_file 4 dfile.rall.1D'[0]' -stim_base 4 -stim_label 6 roll     \
    -stim_file 5 dfile.rall.1D'[1]' -stim_base 5 -stim_label 7 pitch    \
    -stim_file 6 dfile.rall.1D'[2]' -stim_base 6 -stim_label 8 yaw      \
    -stim_file 7 dfile.rall.1D'[3]' -stim_base 7 -stim_label 9 dS       \
    -stim_file 8 dfile.rall.1D'[4]' -stim_base 8 -stim_label 10 dL    \
    -stim_file 9 dfile.rall.1D'[5]' -stim_base 9 -stim_label 11 dP    \
    -stim_file 10 Seed_ts${sd}.1D -stim_label 10 Seed    \
    -stim_file 11 Inter_tsA${sd}.1D -stim_label 11 PPIA    \
    -stim_file 12 Inter_tsB${sd}.1D -stim_label 12 PPIB    \
    -stim_file 13 Inter_tsC${sd}.1D -stim_label 13 PPIC    \
    -rout -tout                  \
    -bucket PPIstat${sd}