Tuesday, August 13, 2013

Psychophysiological Interactions Explained

That last tutorial on psychophysiological interactions - all thirteen minutes of it - may have been a little too much to digest in one sitting. However, you also don't understand what it's like to be in the zone, with the AFNI commands darting from my fingertips and the FMRI information tumbling uncontrollably from my lips; sometimes I get so excited that I accidentally spit a little when I talk, mispronounce words like "particularly," and am unable to string two coherent thoughts together. It's what my dating coach calls the vibe.

With that in mind, I'd like to provide some more supporting commentary. Think of it as supplementary material - that everybody reads!

Psychophysiological interactions (PPIs) are correlations that change depending on a given condition, just as in statistics. (The concept is identical, actually.) For example: Will a drug affect your nervous system differently from a placebo depending on whether you have been drinking or not? Does the effectiveness of Old Spice Nocturnal Creatures Scent compared to Dr. Will Brown's Nutty Time Special Attar (tm) depend on whether your haberdasher is Saks Fifth Avenue, or flophouse? Does the transcendental experience of devouring jar after jar of Nutella depend on whether you are alone, or noshing on it with friends? (Answer: As long as there's enough to go around, let the good times roll, baby!)

So much for the interaction part, but what about this psychophysiological thing? The "psychophysiological" term is composed of two parts, namely - and I'm not trying to insult your intelligence here - a psychological component and a physiological component. The "psychological" part of a PPI refers to which condition your participant was exposed to; perhaps the condition where she was instructed to imagine a warm, fun, exciting environment, such as reading Andy's Brain Blog with all of her sorority sisters during a pajama party on a Saturday night and trying to figure out what, exactly, makes that guy tick. Every time she imagines this, we code it with a 1; every time she imagines a contrasting condition, such as reading some other blog, we code that with a -1. And everything else gets coded as a 0.

Here are pictures to make this more concrete:

Figure 1: Sorority girls reading Andy's Brain Blog (aka Condition 1)
Figure 2: AvsBcoding with 1's, -1's, and 0's
The physiological component, on the other hand, simply refers to the underlying neural activity. Remember, however, that what we see in a typical FMRI timeseries isn't the actual neural activity; it's neural activity that has been smeared (i.e., convolved) with the hemodynamic response function, or HRF. To convert this back to the neural timeseries, we de-smear (i.e., deconvolve) the timeseries by removing the HRF - similar to a Fourier analysis if you have done one of those. (You have? NERD!)

To summarize: we assume that the observed timeseries is a convolution of the underlying neural activity with the HRF. Using NA for neural activity, HRF for the gamma response function, and (x) for the Kroenecker product, symbolically this looks like:

TimeSeries = NA (x) HRF

To isolate the NA part of the equation, we extract the TimeSeries from a contrast or some anatomically defined region, either using the AFNI Clusterize GUI and selecting the concatenated functional runs as an Auxiliary TimeSeries, or through a command such as 3dmaskave or 3dmaskdump. (In either case, make sure that the text file is one column, not one row; if it is one row, use 1dtranspose to fix it.) Let's call this ideal time series Seed_ts.1D.

Next, we generate an HRF using the waver command and shuffle it into a text file:

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

Check this step using 1dplot:

1dplot GammaHR.1D

And you should see something like this:

Now that we have both our TimeSeries and our HRF, we can calculate NA, the underlying neural activity, using 3dTfitter (the -Faltung option, by the way, is German for "deconvolution"):

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

This will tease apart GammaHR.1, from the TimeSeries Seed_ts.1D, and store it in a file called Seed_Neur.1D. This should look similar to the ideal time course we generated in the first step, except now it is in "Neural time":

1dplot Seed_Neur.1D


Ausgezeichnet! All we have to do now is multiply this by the AvsBcoding we made before, and we will have our psychophysiological interaction!

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

This will generate a file, Inter_neu.1D, that is our interaction at the neural level:

1dplot Inter_neu.1D

And, finally, we take it back into "BOLD time" by convolving this timeseries with the canonical HRF:

waver -dt 2 -GAM -input Inter_neu.1D -numout [numTRs] > Inter_ts.1D

1dplot Inter_ts.1D

Let's pause and review what we've done so far:
  1. First, we extracted a timeseries from a voxel or group of voxels somewhere in the brain, either defined by a contrast or an anatomical region (there are several other ways to do this, of course; these are just the two most popular);
  2. We then created a list of 1's, -1's, and 0's for our condition, contrasting condition, and all the other stuff, with one number for each TR;
  3. Next, we created a Gamma basis function using the waver command, simulating the canonical HRF;
  4. We used 3dTfitter to decouple the HRF from the timeseries, moving from "BOLD time" to "Neural time";
  5. We created our psychophysiological interaction by multiplying our AvsBcoding by the Neural TimeSeries; and
  6. Convolved this interaction with the HRF in order to move back to "BOLD time".

Our last step is to enter this interaction and the region's timeseries into our model. Since we are entering this into the model directly as is, and since we are not convolving it with anything, we will use the -stim_file option. Also remember to not use the -stim_base option, which is usually paired with -stim_file; we are still estimating parameters for this regressor, not just including it in the baseline model.

Below, I have slightly modified the 3dDeconvolve script used in the AFNI_data6 directory. The only two lines that have been changed are the stim_file options for regressors 9 and 10:

 #!/bin/tcsh

set subj = "FT"

# run the regression analysis
3dDeconvolve -input pb04.$subj.r*.scale+orig.HEAD                       \
    -polort 3                                                           \
    -num_stimts 10                                                       \
    -stim_times 1 stimuli/AV1_vis.txt 'BLOCK(20,1)'                     \
    -stim_label 1 Vrel                                                  \
    -stim_times 2 stimuli/AV2_aud.txt 'BLOCK(20,1)'                     \
    -stim_label 2 Arel                                                  \
    -stim_file 3 motion_demean.1D'[0]' -stim_base 3 -stim_label 3 roll  \
    -stim_file 4 motion_demean.1D'[1]' -stim_base 4 -stim_label 4 pitch \
    -stim_file 5 motion_demean.1D'[2]' -stim_base 5 -stim_label 5 yaw   \
    -stim_file 6 motion_demean.1D'[3]' -stim_base 6 -stim_label 6 dS    \
    -stim_file 7 motion_demean.1D'[4]' -stim_base 7 -stim_label 7 dL    \
    -stim_file 8 motion_demean.1D'[5]' -stim_base 8 -stim_label 8 dP    \
    -stim_file 9 Seed_ts.1D -stim_label 9 Seed_ts    \
    -stim_file 10 Inter_ts.1D -stim_label 10 Inter_ts    \

    -gltsym 'SYM: Vrel -Arel' -x1D X.xmat.1D -tout -xjpeg X.jpg                             \
    -bucket stats_Inter.$subj -rout



This will generate a beta weights for the Interaction term, which can then be taken to a group-level analysis using your command of choice. You can also do this with the R-values, but be aware that 3dDeconvolve will generate R-squared values; you will need to take the square root of these and then apply a Fisher's R-to-Z transform before you do a group analysis. Because of the extra steps involved, and because interpreting beta coefficients is much more straightforward, I recommend sticking with the beta weights.

Unless you're a real man like my dating coach, who recommends alpha weights.

11 comments:

  1. Just wondering if you need to correct the Seed_Neur.1D file to a minimum of zero before multiplying it by the AvsBcoding.1D file to form the interaction?

    This would avoid the issue of multiplying a negative "Neutral Time" value by a negative condition descriptor (-1) and thus generating a positive value in the interaction file, when it should in fact be a negative value.

    Great blog, by the way!

    ReplyDelete
    Replies
    1. Hi Anonymous,

      The interaction term is looking for voxels where activity is dependent on the condition, which requires weighting positive and negative values for the condition differently; think of a simple [1 -1 -1 1] contrast vector looking for a nonlinear interaction between levels of two conditions. Eliminating all negative values in the neural time file would make it more similar to a linear contrast between the conditions.


      Hope this helps, but let me know if you have any other questions!

      -Andy

      Delete
    2. Hi Andy,

      Thanks for the response!

      I don't think I was clear in my original post. By correcting the Seed_Neur.1D file to a minimum of zero, I meant shifting all values within the file by a certain amount i.e. identifying the most negative value in your file and adding its absolute value to all values in the file, so that the minimum would then be zero.

      The reason I bring this up is because in Karl Friston's original paper on ppi analysis, he mentions correcting the physiological component to a minimum of zero prior to creating the interaction term (p. 223 just above figure 3 - " ... the activity in V1 (corrected to a minimum of zero) ...")*.

      I had presumed this step was taken to avoid the situation I outlined in my original post, but I'm confused as to why this step isn't mentioned in AFNI's guide to ppi.

      Any clarification is most appreciated! :)

      Thanks again!

      *FRISTON, K. J., BUECHEL, C., FINK, G. R., MORRIS, J., ROLLS, E. & DOLAN, R. J. 1997. Psychophysiological and modulatory interactions in neuroimaging. Neuroimage, 6, 218-29.

      Delete
    3. I looked over the paper and talked about it with Gang Chen, who wrote the AFNI guide for doing PPIs, and we were unable to determine why the values were corrected to a baseline of zero. It was unclear from the paper why they did this; your explanation is a possible reason, but again that would be more like a linear contrast than an interaction (because all of condition B would be negative, and all of condition A would be positive, similar to a [1 -1] contrast vector). I may need some clarification from the SPM programmers for why they approach it that way.

      Best,

      -Andy

      Delete
  2. Hi Andy,
    The form of the HRF generated by "waver -dt TR -GAM -inline 1@1 > GammaHR.1D " is like a block response.I think it is because the input Rmodel of your following regression analysis is 'BLOCK(20,1)'.
    I will do PPI analysis to my fMRI data with event-related design. And I think it will be better to choose a response function which is the same as the one used in the following GLM regression. How to generate a HRF for event-related fMRI if the Rmodel used is 'BLOCK5(1,1)' in the following GLM regression?
    I will be appreciated with your answers!
    Thank you!
    Po

    ReplyDelete
    Replies
    1. Hi Po,

      You're right; the waveform should reflect the basis function used in the model. This tutorial was made a while ago, and that was an oversight.

      If you want to create a waveform representing a 5-second boxcar function, then you could do so with 3dDeconvolve:

      3dDeconvolve -polort -1 -nodata 20 1.0 -num_stimts 1 -stim_times 1 '1D: 0' 'BLOCK(5,1)' -x1D Ideal -x1D_stop

      I'm using a 20-second window, but you may want to shorten it or lengthen it as you see fit.

      Best,

      -Andy

      Delete
    2. Hi Andy,
      Thank you for your prompt reply!
      But I think the boxcar function can't match my fMRI data with an event-related design. The Rmodel 'BLOCK5(1,1)' I used in the regression is a basic function for an event-related fMRI design. Maybe the HRF function for my PPI can be generated in this way:
      3dDeconvolve -polort -1 -nodata timepionts TR -num_stimts 1 -stim_times 1 '1D: 0' 'BLOCK5(1,1)' -x1D Ideal -x1D_stop
      But I don't know what's the number of time points should be inputted to the option '-nodata'. Could you tell me how to decide the inputted number of option '-nodata' to generated a HRF can matches the event-related design.
      Po

      Delete
    3. Hey Po,

      You'll want to capture the shape of the entire HRF, which usually goes back to baseline after about 12 seconds. So, assuming you have a TR of 2s, you'd want 7 time points to capture the entire window (2*7 = 14 seconds).

      I notice that you're using the BLOCK5 basis function, with a 1 second duration; is that to reflect the actual duration of the events in your experiment (i.e., they lasted for 1 second)? If so, check your output over a longer time window - say, 30 seconds - and see how long it takes for the HRF of a 1-second event to go back to baseline, and use that time window to generate your HRF.

      Best,

      -Andy

      Delete
  3. Hi Andy,

    first of all thank you very much for your great blog and YouToube channel - they are very helpful for me as an fMRI newbie!!
    I have a question regarding extracting the time series of the seed region in the PPI. I have already selected my seed region I am interested in and, now, I want to extract the time series of that region. But I am quite unsure what GLM I have to use. I usually use a very complex GLM for my first level analysis with lots of different conditions. I used this GLM to get the contrast of interest, i.e. my seed region. Do I have to use this complex GLM again to extract the time series or do I have to create a new GLM including just the conditions I am interested in? (I hope, it becomes clear what my problem is:)

    Thank you!

    ReplyDelete
    Replies
    1. Hello,

      You won't extract the time course from the GLM; the time course is extracted from the preprocessed data that goes into the GLM. For example, if smoothing was your last preprocessing step, then extract the time course from the smoothed data. Then enter that time course as a regressor into your GLM along with your interaction term and the other conditions from your "complex" GLM.

      Best,

      -Andy

      Delete