Generating an Activation Map

Now that the data have been preprocessed, the data are ready for statistical analysis. We will use the General Linear Model, which is  a parametric method to assess the data at each voxel. Let us say that during an experiment we measure a response variable Yj, with j = 1 to J, indexes of observation, and for each observation we have  a set of K( K < J ) explanatory variables, xjk, with k = 1 to K. The explanatory variables may be continuous covariates, functions of covariates, or they may be dummy variables  indicating  the levels of an experimental factor. The general linear model explains the variation in Yj in terms of a linear combination of the explanatory variables, plus an error term.
Yj = xj1b1 + ... + xjkbk+...+ xjKbK + ej

where, bk are unknown parameters, corresponding to each of the K explanatory variables and ej are the errors which are correlated. A description of the procedure is written up as an abstract. After estimating the parameters, the t maps can be calculated using fmrilm. Let us consider the example that was discussed before, which is a block design of a pain study. The analysis is done by starting MATLAB. If you are not very familiar with MATLAB you can take a look at  MATLAB Help Desk , Getting Started, or you can just run demo, in MATLAB. Each functional scan will be analyzed separately.

Before April 2000 we used a different version of fmristat,multistat and tstat_threshold. For a description of the previous version of the software, please visit the old webpage.

The difference between the two versions is that fmristat was completly modified to provide better control over the analysis of the data, and multistat runs faster in certain cases. The basic procedures remain the same, you don't have to reanalyze your data. The results will be a bit different of course, because the assumptions we make are different. The new features are:

In order to analyze the data, fmridesign should be used first. fmridesign is a MATLAB function, so in order to run it , MATLAB has to be opened.

The inputs to this function are the frametime, slicetime, events, S, exclude and the hrf_parameters. The output of the function, X_cache, is a matrix, which will represent the design matrix.

Before running fmridesign, all the variables must be defined:

frametimes is a row vector in which you enter the time, in seconds, at which each frame was acquired. To define this vector, the measurement interval, or the so called "TR" has to be known.The repetition time is set when you do your experiment in the scanner. Be careful, the frame times must be equally spaced (this might be a problem for certain designs of event related  experiments). If for example in the experiment that we are talking about the "TR" was 3s, the frame_times array will look like this:

frametimes=(0:119)*3; or

frametimes =  [0     3    6     9    12    15    18    21   24    27    30    33 ... 351   354  357];

slicetimes is a row vector of relative slice acquisition times. The absolute acquisition time of a slice is frametimes + slicetimes. For different protocols the slice time is different. It also depends on the order in which you do the excitation, ascending, descending or interleaved. If the order is ascending, then the first slice that will be excited will be slice number 1, followed by slice number 2 and so on. If descending, the first slice excited will be your last slice, let's say slice N, the next one will be N-1 and so on. If interleaved, the order is slice 1, slice 3, slice 5, all the odd ones and then slice 2, slice 4, ... all the even ones. You can also pick your own order of excitation, but then you will have to be really careful at the way you calculate your times. With certain protocols, this might be a problem also when you transfer the data from the scanner to the disk on bottom. The header will be incorrect and you will have to rearrange the slices in the correct order. For the protocols we are running in the MNI, the times are:

Mosaic 64;  the excitation order for the Mosaic64 sequence is by default, ascending.

     ~100ms between consecutive slices
    ~140ms from the beginning of the frame to the first slice

BOLD_trig64 (64 x 64 non-mosaic); the excitation order is by default, interleaved.

    is the same as the Mosaic one
     ~100ms between consecutive slices
    ~140ms from the beginning of the frame to the first slice

fMRI_trig (128 x 128 with trigger);  excitation order: ascending (default).

     ~120ms-130ms between slices
     ~140ms-150ms from the beginning of the frame to the first slice

BOLD_fMRI_10sl (128 x 128 no trigger); excitation order : interleaved (default).

     ~120ms-130ms between slices
    The time from slice 1 in frame N to slice 1 in frame N + 1 should be your repetition time (the so called "TR"),  plus  a few  milliseconds.

This means that the acquisition of the slice doesn't start right away when you start your frame, specially for the protocols that wait for a trigger, like the Mosaic 64, the BOLD_trig64 and the fMRI_trig. The first slice is acquired after  a certain time, and this is important , especially for the event-related experiments.
So, for example, in case you are using the 64 x 64, ascending, your slicetime array should look like this, for let's say 25 slices:

slicetimes = [ 0.140  0.240 0.340  0.440 ..... 2.540 ];

If you consider that your experiment is designed in such a way that the delay between the beginning of the frame and the acquisition of the first slice is not important, you can  create an array, slicetimes = ( 0 : ( numslices - 1) ) * 0.1, or for 25 slices let's say:

slicetimes = [0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000  .... 2.3000    2.4000 ];

If you are using the fMRI_trig, the 128 x 128 with the trigger, the slices are excited in ascending order. If you choose to acquire in interleaved order, the first slice would be number 1, then slice number 3, up to slice 25, then slice 2, followed by slice 4, up to slice number 24.
Here is an example: the first slice is excited after 140ms from the beginning of the frame, so the slicetime vector will look like this:

slicetimes = [0.1400    1.7000    0.2600    1.8200    0.3800    1.9400    0.5000 2.0600    0.6200    2.1800    0.7400    2.3000    0.8600    2.4200 0.9800    2.5400    1.1000    2.6600    1.2200    2.7800    1.3400 2.9000    1.4600    3.0200    1.5800];

N.B. In the case you are not synchronizing your stimuli/tasks to the scanner and/or in the case you need precise slice timing information you should consider keeping a log file containing the time when each frame was acquired. The timing information contained in the mincheader is not accurate and the time provided by the scanner is not exactly correct either. Depending upon various acquisition parameters there can be a difference between the requested time between frames and the real time which can be as large as 3%-4%. The best current solution to this is to write your stimulus presentation program in such a way that it records the time when each frame starts. Every time the frame starts the scanner gives a signal (0V to -5V ) which we convert into a middle mouse button click (Logitech PS2 mouse).

events is a matrix whose rows are events and whose columns are:

  1. eventid - an integer from 1 to number of events to identify event type;
  2. eventimes - start of event, synchronized with frame and slice times;
  3. durations- (optional - default is 0) - duration of event;
  4. heights- (optional - default is 1) - height of response for event.
For each event type, the response is a box function starting at the event times, with the specified durations and heights, convolved with the hemodynamic response function (see below). If the duration is zero, the response is the hemodynamic response function whose integral is the specified height - useful for `instantaneous' stimuli such as visual stimuli. The response is then subsampled at the appropriate frame and slice times to create a design matrix for each slice, whose columns correspond to the event id number. events = [] will ignore the events matrix and just use the stimulus design matrix S (see next). Default is [1 0].
Let us now consider our case, which was a block design of three frames rest, three frames hot, three frames rest and three frames warm, all repeated 10 times.The hot event is identified by 1 and the warm event by 2. The event start at times 9,27,45,63... and each has a duration of 9 seconds, with equal height. These will mean we can write: eventid = kron (ones(10,1), [1;2]); eventimes = (0:19)' * 18 + 9; duration = ones(20,1) * 9; height = ones(20,1); events = [eventid eventimes duration height]; or:

events =[

The events can also be supplied by a stimulus design matrix, S. The design matrix consists of rows which represent the frames and the columns are the explanatory (independent) variables of interest, the events. Events are created for each column, beginning at the frame time for each row of S, with a duration equal to the time to the next frame and a height equal to the value of S for that row and column. A constant term is not usually required since it is removed by the polynomial trend terms provided n_poly > 0. The experiment had 120 time frames and 2 conditions hot and warm. For this example, the design matrix consists in a matrix with 120 rows and 2 columns, first one for the hot condition, the second one for the warm condition. At each line, which represents the frames, a 1 or a 0 will be filled in the matrix, depending if the stimulus was on or off. In this case:

S = kron ( ones ( 10,1 ), kron( [0 1 0 0; 0 0 0 1]' , ones ( 3, 1 )))  or :

S = [

N.B. kron stands for Kronecker tensor product, and is a function used in MATLAB. It's used here to give you an example of an easy way to make your design matrices. There are other ways to write the design matrix , just choose the one that is best for you. The one above is just an example. You can make for example one vector called so = ones(3,1) another vector called sz = zeros(3,1) and after that you can combine them and write your design matrix,
S = [sz sz; so sz; sz sz; sz so......].

The two columns of the S matrix represent the two stimuli, the first column is for the hot one and the second column is for the warm stimulus. Each row of the matrix, represents one frame.
One thing to remember is that all values for all frames must be given, because smoothing and lagging  by the heamodynamic function is done BEFORE excluding time points by exclude.

N.B. It is suggested that you use the events matrix to define your events and not the S matrix, even if your experiment is a block design experiment. The results should be the same.

exclude is a list of frames that should be excluded from the analysis. This must be used with Siemens EPI scans to remove the first few frames, which do not represent steady-state images. Default is [], but I would recommend that you exclude the first 3 scans.

exclude = [1 2 3];

hrf_parameters is a matrix whose rows are six parameters for the hemodynamic response function; one row for each event type and column of S (if there is just one row, this is repeated as necessary). The hrf is modeled as the difference of two gamma density functions: "Deconvolution of impulse response in event-related BOLD fMRI." (Glover, NeuroImage, 9:416-429). The components of HRF_PARAMETERS are:

  1. PEAK1: time to the peak of the first gamma density;
  2. FWHM1: approximate FWHM of the first gamma density;
  3. PEAK2: time to the peak of the second gamma density;
  4. FWHM2: approximate FWHM of the second gamma density;
  5. DIP: coefficient of the second gamma density;

    Final hrf is: gamma1 / max(gamma1) - DIP * gamma2 / max(gamma2) scaled so that its total integral is 1.

  6. FIT_SCALE: 1 - fit the time scale of the hrf by convolving its derivative with the specified column of the design matrix, to create an additional column for the design matrix. Dividing the effect of this column by the effect of the hrf itself estimates the scale shift. 0 ignores this option.
If PEAK1 = 0 then there is no smoothing of that event type with the hrf. Default is: [5.4 5.2 10.8 7.35 0.35 0] chosen by Glover (1999) for an auditory stimulus.

For the example we are considering, we will choose the default values:

hrf_parameters = [5.4 5.2 10.8 7.35 0.35 0];

You could take a look at the hemodynamic response function, just by doing this:

frametimes= (0 : 119) * 3;

X = fmridesign(frametimes, 0, [1 0], [ ], [ ], hrf_parameters);

plot(frametimes, X)

Now, that all the variables are defined,the final design matrix can be calculated, just by calling the fmridesign function in MATLAB:

X_cache = fmridesign (frametimes, slicetimes, events, S, exclude, hrf_parameters)

The output of this function, the matrix X_cache will be used to calculate the stats for the experiment we designed. For that you will have to use the function fmrilm. fmrilm reads in a fMRI file and writes out a T statistic file. The inputs to fmrilm are: input_file, output_file_base, X_cache, contrast, exclude, which_stats, fwhm_rho, n_poly. The output, is a matrix with two elements: DF, the residual degrees of freedom of the statistics and P, the degrees of freedom of the set of contrasts

Again the variables have to be defined first, before we can use the function.

input_file, is the minc file we would like to analyze and it is the dependent variable. In our case:

input_file = 'smith_john_19971024_1_105235_mri_MC.mnc'

output_file_base, is a matrix whose rows are the base for output statistics, one base for each row of contrast, padded with extra blanks if necessary; they will be removed before use.
Let us consider that in our example we would like to compare hot minus warm (hmw) and an overall combined effect of hot plus warm (hpw) stimulus versus the baseline. So the output_file_base will look like this:

output_file_base =['smith_john_19971024_1_105235_mri_MC_hmw'; 'smith_john_19971024_1_105235_mri_MC_hpw'];

X_cache is the design matrix for all slices and it is usually the output from fmridesign. Rows are the non-excluded frames, columns are all the regressor variables for slice 1, then slice 2, i.e. with slices running slowest. If X_cache only contains variables for just one slice, it is applied to all slices. You can also provide your own design matrix, if you like, and in this case you will not have to run fmridesign. For our example, we will consider the X_cache, that we just calculated before.

contrast is a matrix whose rows are contrasts for the statistic images, with row length equal to the number regressor variables in X_CACHE plus n_poly+1 ( n_poly will be described later ); the extra values allow you to set contrasts in the drift terms as well as the regressor variables. As we disscussed, we want to test for hot minus warm and also for hot plus warm, against the baseline, for this case, the contrast matrix will look like:

contrast = [1-10000 ;
110000] ;

exclude is a list of frames that should be excluded from the analysis. This must be used with Siemens EPI scans to remove the first few frames, which do not represent steady-state images. Default is [1], but we suggest that you take out the first three images.

exclude = [1 2 3];

which_stats is a logical matrix indicating output statistics by 1. Rows correspond to rows of contrast, columns correspond to:

  1. T statistic image, output_file_base_tstat.mnc. The degrees of freedom is DF. Note that tstat=effect/sdeffect.
  2. effect image, output_file_base_effect.mnc
  3. standard deviation of the effect, output_file_base_sdeffect.mnc.
  4. F-statistic image, output_file_base_Fstat.mnc, for testing all rows of contrast simultaneously. The degrees of freedom are [DF, P].
  5. the AR parameter, output_file_base_rho.mnc.
  6. the residuals from the model, output_file_base_resid.mnc, only for the non-excluded frames (warning: uses lots of memory).
  7. the whitened residuals from the model, output_file_base_wresid.mnc, only for the non-excluded frames (warning: uses lots of memory).
If which_stats is a row vector of length 7, then it is used for all contrasts; if the number of columns is less than 7 it is padded with zeros. Default is 1. In our example:

which_stats =[1111111];

fwhm_rho is the fwhm in mm of a 3D Gaussian kernel used to smooth the autocorrelation. This is NOT the fwhm of the image data as used by tstat_threshold, but should be fairly large. Setting it to Inf smooths the autocorrelation to 0, i.e. it assumes the scans are uncorrelated (useful for TR>10 seconds). Default is 15.

fwhm_rho = 15;

n_poly is the order of polynomial trend terms to remove.
If n_poly > 0, it will remove a polynom with the order of n_poly.
If n_poly = 0, it will remove just the constatnt term and no trends.
If n_poly < 0, it will not remove anything, in which case the design matrix is completely determined by eventimes and S.
The  default value for n_poly is 3. This will remove the scanner drift.

n_poly = 3;

At this point all the variables are set, and the command can be issued in MATLAB:

[DF, P] = fmrilm( input_file, output_file_base, X_cache, contrast, exclude, which_stats, fwhm_rho, n_poly);

After running fmrilm, the statistics of the image is calculated, and the *tstat.mnc file is the t map file. The statistical significance thresholds for the t-maps should be calculated. To do that tstat_threshold should be used, also in MATLAB.
Previous Next Back to table of contexts