df = fmristat ( INPUT_FILE , FRAME_TIMES , SLICE_TIMES , EVENT_TIMES , S, CONTRAST , EXCLUDE , OUTPUT_FILE_BASE , NUM_STATS , MEAN_LAG , SD_LAG , FWHM_RHO , N_POLY )
INPUT_FILE is the minc file we would like to analyze and it is the dependent variable.
input_file = 'smith_john_19971024_1_105235_mri_MC.mnc'
N.B. I will use capital letters for defining the variables. When you use them in MATLAB, you can use either lower case or capital case, as long as you are consistent.
FRAME_TIMES is a row vector
in which you enter the time at which each frame was acquired. If a single
scalar is supplied, this is used as the frame time separation, starting at
zero. This value will be the measurement interval, or the so called "TR",
that you 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];
SLICE_TIMES is a row vector
of relative slice acquisition times. If a single scalar is supplied, this is
used as the slice time separation, starting at zero. For different protocols
the slice time is different. It also depends on the order 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 really be careful at the way
you calculate your times. 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: interleaved (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
enter just slicetimes = 0.1 and that will 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 excitied in
interleaved order, starting with slice number 1, then slice number 3, up to
slice 25, then slice 2, followed by slice 4, up to slice number
24.
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];
EVENT_TIMES is a matrix of event times this is to be used when an event related experiment is designed ( an example will be given later ). For the block design experiment , EVENT_TIMES should be EVENT_TIMES=[] , this will ignore the EVENT_TIMES and the stimulus design matrix, S will be used.
The stimulus design matrix, S, has to be written. The design matrix consists of rows which represent the frames and the columns are the explanatory (independent) variables of interest. 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.
If your experiment is event related, then you just have to make S, an empty
matrix, S = [], this will ignore the stimulus design matrix and will
consider only the EVENT_TIMES.
CONTRAST is a row vector of contrasts for the T statistic image with a length equal to the number of columns of [ EVENT_TIMES S ]. Default is [ 1 0 ... 0 ], i.e. it picks out the first column. For our particular example, contrast is :
contrast = [1 -1]; that means it will tell you which linear combination of the columns of S you are interested in. Here it is hot-warm, the first column, hot, minus the second column, warm.
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], i.e. the first frame is excluded, but I would recommend that you exclude the first 3 scans.
exclude = [1 2 3];
OUTPUT_FILE_BASE: Base for output statistics. Default is the first INPUT_FILES name minus the .mnc extension, or minus the .mnc.gz extension if it is compressed. For the example mentioned before
output_file_base =
'smith_john_19971024_1_105235_mri_MC', or you can choose any name you
want as long as you don't overwrite files and you know which t-map
corresponds to a certain file.
To get the stats: tstats, effects and their sd's, you will have to specify the NUM_STATS :
NUM_STATS: output the
following images up to NUM_STATS, the number you input:
1: T statistic image,
OUTPUT_FILE_BASE_tstat.mnc. The degrees of freedom is fmristat_df, which is
printed out at the end.
2: effect image,
OUTPUT_FILE_BASE_effect.mnc
3: standard deviation of the
effect, OUTPUT_FILE_BASE_sdeffect.mnc.
4: the AR parameter,
OUTPUT_FILE_BASE_rho.mnc.
NB:
tstat=effect/sdeffect.
Default value for NUM_STATS is
1.
In our case, let's say we would like to be able to calculate the stats for a group of subjects, and also to calculate the stats for the same subject, but different sessions. This is why we need the tstats, the effects and their sd's, thus num_stats = 3. If you would like to have the autocorrelation map, then:
num_stats = 4;
MEAN_LAG and SD_LAG are parameters for the hemodynamic response function. This is modeled as a gamma density function with mean lag equal to MEAN_LAG, and standard deviation equal to SD_LAG. Defaluts are MEAN_LAG=6, SD_LAG = 3 (seconds), taken from Zarahn et al. (1997). Setting MEAN_LAG = 0 will remove all lagging. If your experiment if designed in such a way that you don't need to take in account the heamodynamic response function, then you should set MEAN_LAG=0.
The peak lag is less than the MEAN_LAG:
peak lag = MEAN_LAG - SD_LAG2 / MEAN_LAG.
If you do not enter any
value the default value will be used.
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. Default is 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 = fmristat(input_file,frametimes, slicetimes, [], S, contrast, exclude, output_file_base,num_stats);
in the case where you kept the default values, if you changed some of the variables then you will have to enter the command in MATLAB with the new values:
df = fmristat ( input_file, frametimes , slicetimes , [] , S, contrast , exclude , output_file_base , num_stats , mean_lag , sd_lag , fwhm_rho , n_poly )
After running fmristat you
will find four new files in your working directory:
smith_john_19971024_1_105235_mri_MC_effect.mnc,
smith_john_19971024_1_105235_mri_MC_sdeffect.mnc,
smith_john_19971024_1_105235_mri_MC_tstat.mnc,
smith_john_19971024_1_105235_mri_MC_rho.mnc.
Besides the four files, the number of degrees of freedom will be printed at the end of the fmristat program. You should remember this number, it will be used to estimate the threshold for the t-maps.
Till now, the analysis of a block design experiment was discussed, but
an event related experiment can be
designed also. Consider that you designed your experiment to have 120
frames, equally spaced, 3s "TR" , and the hot stimulus was given once at the
beginning of scan 4, 16, 28, ... etc. and that the warm stimulus was given
at scans 10, 22, 34, ...
The only difference between the block design
and the event related will be that you will have to provide times of the
events for each type of stimulus, instead of a design matrix S. The rest is
the same. If the frames take place at the times shown in the next table,
then the eventimes matrix should be:
| Frame # | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | ... |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Times (s) | 0 | 3 | 6 | 9 | 12 | 15 | 18 | 21 | 24 | 27 | 30 | 33 | 36 | 39 | 42 | 45 | 48 | 51 | 54 | ... |
eventimes = [
9 27
45 63
81 99
117 135
153 171
189 207
225 243
261 279
297 315
333 351];
If in your experiment you have more events of one type than another, pad the matrix with any time beyond the last frame time, e.g. 500.
After you provided the event times matrix, the command in Matlab will be:
df = fmristat(input_file, frametimes, slicetimes, eventimes, [], contrast, exclude, output_file_base, num_stats);
NB: You can combine both eventimes and S if you like.
As it was mentioned before the t maps are
obtained by running fmristat.m. 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.m should be used, also in
Matlab.
| Previous | Next | Back to table of contexts |
|---|