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 ,

Before April 2000 we used a
different version of ** fmristat**,

The difference between the
two versions is that ** fmristat** was completly modified to provide better
control over the analysis of the data, and

was split into two:*fmristat*, which will set up the design matrix.*fmridesign*, which will do the fitting*fmrilm*

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

The inputs to this function are the * frametime*,

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

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:

**eventid**- an integer from 1 to number of events to identify event type;**eventimes**- start of event, synchronized with frame and slice times;**durations**- (optional - default is 0) - duration of event;**heights**- (optional - default is 1) - height of response for event.

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

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

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

* 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:

- PEAK1: time to the peak of the first gamma density;
- FWHM1: approximate FWHM of the first gamma density;
- PEAK2: time to the peak of the second gamma density;
- FWHM2: approximate FWHM of the second gamma density;
- 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.

- 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:

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

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

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'];

contrast = [ | 1 | -1 | 0 | 0 | 0 | 0 ; | |

1 | 1 | 0 | 0 | 0 | 0 | ] ; |

exclude = [1 2 3];

- T statistic image,
**output_file_base_tstat.mnc**. The degrees of freedom is DF. Note that tstat=effect/sdeffect. - effect image,
**output_file_base_effect.mnc** - standard deviation of the effect,
**output_file_base_sdeffect.mnc**. - F-statistic image,
**output_file_base_Fstat.mnc**, for testing all rows ofsimultaneously. The degrees of freedom are [DF, P].*contrast* - the AR parameter,
**output_file_base_rho.mnc**. - the residuals from the model,
**output_file_base_resid.mnc**, only for the non-excluded frames (warning: uses lots of memory). - the whitened residuals from the model,
**output_file_base_wresid.mnc**, only for the non-excluded frames (warning: uses lots of memory).

which_stats = | [ | 1 | 1 | 1 | 1 | 1 | 1 | 1 | ]; |

fwhm_rho = 15;

If

If

If

The default value for

n_poly = 3;

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

Previous |
Next |
Back to table of contexts |
---|