# Context-Dependent Correlation Analysis or Generalized PPI

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^{2}, not correlation coefficient R itself. So we need to take square root of R^{2 }and find out its sign based on the sign of its corresponding beta value:

*3dcalc -a *Corr_subj01+orig'[SubbrickForR^{2}]'* -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.

Last modified 2015-01-07 15:53