Generating an Activation Map with fmristat

Now that the data have been preprocessed, we 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 ) response variables, xjk, with k = 1 to K. The response 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 response 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. Each functional scan will be analyzed separately.

Overview of commands used in fmristat

The 3 basic functions that you will use in fmristat are the following:

The analysis is done in Matlab which is run from the command line by typing in 'matlab' (without the quotes). If you are not very familiar with Matlab you can take a look at Matlab's Help Desk on the web, or you can just run demo, from within Matlab. The basic idea behind the Matlab commands is simply (a) specifying the variables that will be used in each of the three functions above, and (b) executing these functions with the specified variables.

fmridesign

We begin with fmridesign. The inputs to this function are the frametimes, slicetimes, and events. The output of this function, X_cache, is a struct, which is used in fmrilm.

Before running fmridesign, we must define the following variables:

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 "TR" must to be known. The repetition time is set when you do your experiment in the scanner. The TR will normally be constant throughout the run, unless you are using physiological gating at the scanner. Let us say for this example, the TR is 3s; the frametimes variable will look like this (Note: the first frame should be set to t=0s):

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 slicetimes are ordered in the row vector according to their voxel coordinates. For axial slices, the z-space coordinate value of 0 would refer to slice #1 and the time for this slice would be entered into the first column of slicetimes. If your minc image dimension steps are negative, the first slice will be the topmost, superior slice - this is the case for most fmri images acquired at the BIC. Once you have determined where the first slice can be found, the slicetimes variable would follow the convention:

slicetimes = [  <time for slice #1>   <time for slice #2>  ...  <time for slice #N>  ]

The slicetimes will depend on the delay in TR that was specified for your protocol at the scanner. It will also depend on the MR sequence that was used in your scans, and the slice excitation order: ascending, descending or interleaved.

A. Slice Timings

You can determine your slice time interval with the following calculation:

slicetime = [TR - delay in TR]  /  [total number of slices]

NB: There may be a very small delay (12 to 14ms) from the acquisition of the first slice and the trigger delivery from the scanner. This should not impact your analysis. Overall, you can expect the total scan time of a complete run to be accurate to within +/- 10 milliseconds.

B. Slice Excitation Direction

1. Ascending

If the order is ascending, then the first axial slice that will be acquired will be the most inferior slice (toward the the subject's feet), and the last slice will be the most superior slice (toward the top of the head).

If for example, you are using the mosaic 64, ascending, your slices are acquired inferior to superior (the foot to head direction). Therefore the top slice is the last one to be imaged and the slicetimes variable will look like this for 25 slices:

slicetimes = [  2.40  2.30  2.20 ... 0.10  0  ];

2. Descending

If descending, axial or transverse slices will be acquired superior to inferior, or head to foot, with the first slice being near the top of the brain and the last slice being the most inferior one.

The descending mode is the default ordering, but you should confirm this for your protocol with the MR technician. If your experiment is performed with the mosaic 64 sequence descending, your slicetimes might look like the following:

slicetimes = [  0  0.10  0.20 ... 2.30  2.40  ];

You could make use of whatever Matlab shortcuts you may know. You could for example create an array with a numslices variable being the number of slices in your functional images:

slicetimes = ( 0 : ( numslices - 1 ) ) * 0.100;

3. Interleaved

The order of acquisition of an interleaved scan depends whether the total number of slices is odd or even.

A. Odd numbered total amount of slices
If your fMRI images were acquired with an odd numbered total amount of slices, the default interleaved sequence begins with the most inferior slice and acquires every odd numbered slice first, ascending to the top of the brain. The acquisition then returns to the lowest even slice and ascends to the top of the brain acquiring all the even slices.

As an example, if you are using the 128 matrix fMRI sequence, and you choose to acquire in interleaved order, the lowest axial slice is excited first, and the second slice from the bottom of the brain will be acquired only after all the odd numbered slices are acquired. So the slicetimes vector will look like the following for an odd numbered 15 slice acquisition with 110 ms per slice:

slicetimes = [  0.7700   1.5400   0.6600   1.4300   0.5500   1.3200   0.4400 ...
   1.2100   0.3300   1.1000   0.2200   0.9900   0.1100   0.8800   0  ];  

The shortcut Matlab code that generated these interleaved slicetimes is a bit tricky and there are several ways to do it. Here is one way:

slicetimes_tmp = (14:-1:0) * 0.110;
s1 = slicetimes_tmp(8:15);
s2 = slicetimes_tmp(1:7);
slicetimes(1:2:15) = s1;
slicetimes(2:2:15) = s2;

B. Even numbered total amount of slices
If your scanning protocol specifies an even numbered total amount of slices, the default interleaved sequence starts with the even slices acquired first, then the odds. In this case your order of acquisition would be first slices 2, 4, 6, ..., then slices 1, 3, 5, ...

As an example with the same slice parameters as above, except with an even numbered total slices, say 14, we have the following:

slicetimes = [  0.6600   1.4300   0.5500   1.3200   0.4400   1.2100   0.3300 ...
   1.1000   0.2200   0.9900   0.1100   0.8800   0   0.7700  ];  

Note: If you require accurate timing information, you should try to configure your stimulus presentation software to record the time when each trigger signal is received, even if you do not require the trigger to advance the stimulus presentation. You may recall that each time a frame begins, the scanner emits a TTL pulse which we convert into keyboard keypresses or mouse clicks. Refer back to the fmri resources section for a review.


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, or HRF. If the duration is zero, the response is the HRF whose integral is the specified height - this is 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.

Let us now consider our example case, which was a block design of three frames rest, three frames hot, three frames rest and three frames warm, all repeated 5 times. The hot event is identified by 1 and the warm event by 2 and the rest condition is left unspecified since it is the baseline. The events start at times 30, 90, 150, ... and each has a duration of 30 seconds, with equal height. This means we can write:

eventid = kron (ones(3,1), [1;2]);
eventimes = (0:5)' * 60 + 30;
duration = ones(6,1) * 30;
height = ones(6,1);
events = [ eventid eventimes duration height ];

or:

events = [
130301
290301
1150301
2210301
1270301
2330301 ];

Note: "kron" stands for Kronecker tensor product, and it is a function used in Matlab. The use here is just to give you an example of an easy way to make your design matrices. You are free to write the events matrix in any way that you prefer.

One thing to keep in mind is that all values for all frames must be given, because the smoothing and lagging by the hemodynamic function is done on all frames before fmrilm excludes the initial time points.

At this point, all of the necessary variables are defined. There are other parameters which you can specify - type 'help fmridesign' for a complete list. If you choose not to specify any trailing variables, you will obtain their default values. The final design matrix can be calculated by calling the fmridesign function in Matlab and equating its output to an X_cache variable:


X_cache = fmridesign (frametimes, slicetimes, events)

The output of this function, the struct X_cache will be used as an input for the next step in the analysis, the function fmrilm.

fmrilm

This fmristat function reads in an 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_cor, and n_ntrends.

There are several values output, including df.t, the effective degrees of freedom of the T statistic. And if an F-statistic is performed, df.F is also available. Again the input 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 = 'study_subj001_20031019_111740_3_mri_MC.mnc'

output_file_base, is a matrix whose rows are the base filenames for the output statistics, one filename base for each row of contrast is needed. 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 = ['study_subj001_20031019_111740_3_mri_MC_hmw';
 'study_subj001_20031019_111740_3_mri_MC_hpw'];


X_cache is the design matrix for all slices and it is the output from fmridesign.

contrast is a matrix whose rows are contrasts for the response variables, and columns are the contrasts. The column numbers will correspond to the event ids that were used in the events matrix above. Since in this example, we want to test for hot minus warm, and also for hot plus warm, against the baseline, the contrast matrix will look like:
contrast = [ 1-1;
 1 1];

exclude is a list of frames that should be excluded from the analysis. This could be used to remove the first few frames of your functional run, which may not represent steady-state images. You can also use this variable to exclude frames during which your subject moved more than 1mm or 1 degree.

exclude = [  1  2  ];


which_stats is a matrix that indicates desired output statistics with a series of strings. Rows correspond to the rows of contrast. The required strings to use are as follows:
If which_stats is a row vector, then it is used for all contrasts. The default output is _mag_t, which gives you just the t-stat output. Do not forget to include the leading underscore character:

which_stats = '_mag_t  _fwhm  _mag_sd  _mag_ef' ;

fwhm_cor is the fwhm in mm of a 3D Gaussian kernel used to smooth the autocorrelation of residuals. This is NOT the fwhm of the image data as used by stat_threshold, but should be fairly large. Setting it to Inf smooths the autocorrelation to 0, i.e. it assumes the scans are uncorrelated (this is appropriate when the TR>10 seconds). If this value is negative it is taken to be the desired df. The default is -100, to target a df of 100.

fwhm_cor = -100;

n_trends is the number of trends to remove. This is a 3 column vector of the following form:
n_trends=[n_temporal n_spatial pcnt]
The  default value for n_trends is [  3  1  1  ]. This will convert the data to percent of whole volume and spatially and temporally detrend it.

Note: If you wish to execute the analysis according to the method employed by previous versions of fmristat prior to 2002, you should set n_trends = [  3  0  0  ];

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


[df spatial_av] = fmrilm (input_file, output_file_base, X_cache, contrast, exclude, which_stats, fwhm_cor, n_trends);


After running fmrilm, the statistics are calculated, and the output_file_base_mag_t.mnc file is the t-stat map file. The statistical significance thresholds for the t-maps should be calculated. To do that, stat_threshold is used, also on the Matlab command line.



Previous Back to table of contents Next