Generating an Activation Map

where, b_k are unknown parameters, corresponding to each of the K explanatory variables and e_j 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:

• fmristat was split into two:
  1. fmridesign, which will set up the design matrix.
  2. fmrilm, which will do the fitting
• the hemodynamic response function was defined in a different way, it includes a dip after the peak;
• the events are now specified by time, duration and height of response;
• in order to reduce computations, multiple contrasts can be provided;
• F-tests for all the contrasts will be calculated;
• there is a greater control over which statistics should be out putted;
• multistat is faster when no smoothing of the random effect is required, or when the design matrix has the same number of rows and columns.

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.

events =[

	1	9	9	1
	2	27	9	1
	1	45	9	1
	2	63	9	1
	1	81	9	1
	2	99	9	1
	1	117	9	1
	2	135	9	1
	1	153	9	1
	2	171	9	1
	1	189	9	1
	2	207	9	1
	1	225	9	1
	2	243	9	1
	1	261	9	1
	2	279	9	1
	1	297	9	1
	2	315	9	1
	1	333	9	1
	2	351	9	1	];

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 = [

	0	0
	0	0
	0	0
	1	0
	1	0
	1	0
	0	0
	0	0
	0	0
	0	1
	0	1
	0	1
	0	0
	0	0
	0	0
	1	0
	1	0
	1	0
	............
	];

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.

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:

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	-1	0	0	0	0 ;
	1	1	0	0	0	0	] ;

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).

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

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