3dDeconvolve

 

INTRODUCTION

 

3dDeconvolve is a versatile command line program for AFNI.  There is extensive documentation for this program in the file 3dDeconvolve.pdf or 3dDeconvolve.ps.  I found those manuals a bit dry and most examples are theoretical and do not refer to actual experiments.  So, I wrote this little ditty and hopefully after coursing through the contents of this and other related manuals you will have an idea of the tools you need to use this program to analyze some data.

 

Here are a few unique capabilities of 3dDeconvolve.  3dDeconvolve can be used to assess an experimental design before participants are scanned.  In this manner several competing designs can be compared.  The program also allows one to calculate an impulse response function for each voxel.  Most importantly, 3dDeconvolve has the capability of calculating statistics from datasets that have greater than one stimulus.  The different stimulus can then be compared in terms of the area under the IRF curve as well as maximum response.  It is this capability that may prove useful in more complex experimental design.

 

For information on the theory and more detailed descriptions of options see 3dDeconvolve.pdf or 3dDeconvolve.ps (type these into a search engine and lots of different sites have links to these files)

 

I'm going to deal with different options as I go along otherwise this manual will be just as dry as the original.

 

USE 1:  Evaluating experimental design

 

So let's examine the experimental design of a block design containing 4 blocks of 8 stimulus and 8 rest scans.  That will give us a total of 64 volumes.  In order to examine this set-up we need to write a 1D file containing the stimulation file listed above.  This file contains 1 column with 64 rows with ones where the stimulation will be and 0s for the rest phase.  Look at the file Example1a.1D

 

We now run the program using the following code

 

3dDeconvolve              \

-nodata                                    \

-nlast 63                                   \

-polort 0                                   \

-num_stimts 1               \

-stim_file 1 Example1a.1D            \

-stim_maxlag 1 7

 

This prints out the following information (I removed the inverse X matrix..call me if you want to see it)

 

Program:          3dDeconvolve

Author:           B. Douglas Ward

Initial Release:  02 Sept 1998

Latest Revision:  06 July 2001

 

Stimulus: Stim #1 

  h[ 0] norm. std. dev. =   0.5179 

  h[ 1] norm. std. dev. =   0.5345 

  h[ 2] norm. std. dev. =   0.5345 

  h[ 3] norm. std. dev. =   0.5345 

  h[ 4] norm. std. dev. =   0.5345 

  h[ 5] norm. std. dev. =   0.5345 

  h[ 6] norm. std. dev. =   0.5345 

  h[ 7] norm. std. dev. =   0.5179 

 

OK we see then the the norm. std. dev. for any given lag is approximately 0.53.

 

Now what do the different options in the 3dDeconvolve command line mean?

 

-nodata                        his means that 3dDeconvolve will be evaluating an experimental design                                       not measured data

-nlast                            specifies the last image to be used in the analysis...default is the number                                     of images in the dataset

-polort                          the degree of the polynomil baseline model.  0 means there is no linear                                            trend.  0 is used for evaluating experimental design and 1                                                     for a model with linear trend (3d+time datasets)

-num_stimts                  number of input stimulus time series to be used

-stim_file k                   stim_file is the name of the stimulus file and k corresponds to the                                            number of the stimulus.  In the above example k=1 as there is only one                                           stim file.

-stim_maxlag k n          n specifies the maximum timelag to be used when estimating IRF that                                           corresponds to the kth stimulus file

 

For more examples of using 3dDeconvolve to assess experimental design see the document Experimental design.doc.

 

 

USE 2  Estimating the HRF for a given stimulus and calculating statistics for a one stimulus dataset!

 

For this example I will use the dataset C-E-Mad to calculate hrfs and simple statistics.  Now  C-E-Mad is a simple block design with 5 equilibration scans and then 8 blocks with 6 active scans followed by 13 conrol.  The total number of scans is 5 + 19X8 = 157.

 

I wrote out a 1D file of the name mad_sad_glad.1D that contains one column wit 157 rows.  In rows that correspond to an active period the value is 1.  At all other rows the value will be zero.

 

The following was entered in the command line…

 

akasha% 3dDeconvolve \

-input C-E-Mad+orig \

-polort 1 -nfirst 5 \

-num_stimts 1 \

-stim_file 1 mad_sad_glad.1D  -stim_label 1 mad_sad_glad \

-stim_minlag 1 0 -stim_maxlag 1 8 \

-iresp 1 C-E-Mad.irf \

-tout -rout -fout -bucket C-E-Mad_nf5.bucket

 

New Options Introduced

 

-input                           used to specify the 3d+time dataset the model will be matched to

-stim_label 1                 provides a prefix for stimulus number one that will appear in the bucket dataset

-stim_minlag 1 0            this is the default condition so this is sort of redundant

-iresp 1                        produces a dataset that estimates the IRF at each voxel.  This dataset will contain 9 points for each voxel.  TR 0 (minlag), 1, 2, …, 8 (maxlag)

-tout                             T-statisics are outputfor each regression time lag

-rout                             R2 is output for the full and partial model

-fout                             F-stats are output for full and partial model

 

Notice in this example that the partial and full models are the same as there is only one input stimulus.

 

In any event the following will be displayed on the monitor as the computer calculates

 

Program:          3dDeconvolve

Author:           B. Douglas Ward

Initial Release:  02 Sept 1998

Latest Revision:  06 July 2001

 

Writing `bucket' dataset into ./C-E-Mad_nf5.bucket+orig.HEAD

--- Writing 3d+time dataset into ./C-E-Mad.irf+orig.HEAD

 

Run AFNI to see the results of your beautiful work.

 

To view how your model fits the data.  Load in ax as the anatomical and load in C-E-Mad_nf5 as your functional dataset.  This functional dataset is a bucket type dataset with 26 sub-bricks.

 

Quickly here is a quick rundown of what they are

0:  co-efficient for  constant nopise value

1:  t-stat for above

2:  co-efficient for linear trend

3:  t-stat for above

Briks 4+5 are for a coe-fficient and t-stat for the data at timelag 0

Briks 6+7 are for a co-efficient and t-stat for data at timelag 1

This pattern continues for Briks 8-21

The most useful sub-bricks follow 22 is for the R2 statistics corresponding to stim 1 (mad_sad_glad)

Whereas, sub-brick 23 is for parital F-stats fro mad_sad_glad

Sub-bricks 24+25 are identical to the above as they are statistics for the whole model and since we had just one stimulus the whole model is the partial model.

 

HOW TO VIEW THE RESULTS

 

I find that the best way to view the results is to load #22 (#24) as the func dataset.

Then load #23 (F-stat) as the Thr Brik.

 

Move the F-test slider bar to a value of around 6.00 or to any other desired value.  In terms of the intensity slider I find that segmenting the data as follows gives a pretty picture

0.0-0.20 nothing

0.2-0.35 red

0.35-0.5 orange

0.5-1.00 yellow

 

The following is a picture I made with an F-value cut-off of 6.0.

 

C-E-Mad_F6_nf5.pnm

 

So there you are able to make a map that corresponds to the BoldFOLD maps you are used to.

 

 

 

LOOKING AT THE HRF FOR A GIVEN VOXEL:

 

Using the –iresp option in the 3dDeconvolve command line we were able to create a 3d+time dataset that contains the fitted impulse reponse function at every voxel.  The dataset will contain 9 subricks as this was the number of timelags used in the command line…minlag 0 and maxlag 8.  So I guess you want to take a peek and look at the results.

 

To look at this dataset just load it as an anatomical image.  It is best to use an overlay to see which voxels are active during the stimulus.  So load in the irf timecourse dataset.  It will have the suffix .irf.  Now hit the graph button to see the irf at every voxel.

 

I printed out pnms for a couple of voxels that were particularly nice.

 

Take a look at

X66y73z1-nf5.pnm as well as

X67y70z1-nf5.pnm

 

So there you have it.  You can examine the irf and see that it corresponds nicely to the traditional hrf.  There is a small undershoot in these functions that occurs a second or two after the stimulus…likely corresponding to the undershoot.  Notice also the peak of the response is around 2-3TRs so about 2.5 X 1.65 = 4 seconds or there abouts.  Then the function falls off. 

 

An interesting feature here is that you can use the option write center to print out a 1d timeseries for a voxel and you can use this as a reference waveform.  Moreuseful may be using the .irf dataset with a ROI (Brodmann area etc.) and writing out an average time series file for the ROI and then using this in experimental analysis.  Once again just the thoughts and ideas of a lowly summer student so take what you see as useful.