Ý
fMRI Analysis Package
by Pawel Skudlarski

The Official Manual

last update

October 1998

Ý

0. Table of Contents

1. Introductory notes
2. Basic Physics and Physiology of fMRI
3. Structure of this program
4. Typical steps in the analysis of a functional project
5. How to get started on SUN or SGI workstation
6. Setup File - description of the study
7. Running read_fmri
8. Additional data analysis
9. Correlation analysis
10. Fourier analysis
11. Running a Batch Job
12. Viewing Saved Results
13. Registration of images into the Talairach Atlas
14. Combining data from several studies
15. Postprocessing Commands
16. Dictionary of useful commands
17. Using optical disks
18. Some technical details
19. Troubleshooting (before you ask Pawel)
20. short description of some statistics
21. ROI time course

Useful links

Questions and comments and BUUUUGS!!!!!!

Please send e-mail to pawel.skudlarski@yale.edu. Or call him at 785-5462.
Back to the Table of Contents

1. Introductory notes

1.1. OnLine help

Most of the MATLAB and fMRI comands gives a short help entry if you type: help command_name. To find location of the MATLAB file defining this command type: which command_name.

1.2. About this program

This program was written by Pawel Skudlarski with help of Cheryl Lacadie, Todd Constable and Adam Anderson in the NMR Research group at Department of Diagnostic Radiology of Yale Medical School. It is being used by this group as a main tool for processing of the functional brain imaging data. With any questions you may contact Pawel Skudlarski at:

Department of Diagnostic Radiology

pawel.skudlarski@yale.edu
Back to the Table of Contents

2. Basic Physics and Physiology of fMRI

Functional MRI (Magnetic Resonance Imaging) is a method of tracking the brain activation's by using multiple MRI images of the brain. The differences between images taken while subject is performing different tasks tell as about the regions of the brain that activated. There is still a lot to be learn about the link between the brain activation's and changes in the MRI but currently the mostly accepted is the BOLD (Blood Oxygen Level Dependent) mechanism. The local increase of the brain activation increases the metabolism rate and the oxygen consumption. This causes increase of the large blood flow increase that overcompensates for the increased use of oxygen. In result the blood in the active region becomes more oxygenated. This means that in contains more oxyhemoglobin and less deoxyhemoglobin. Luckily for MRI requirement these two forms of hemoglobin vary in its magnetic properties. Deoxyhemoglobin is paramagnetic. The higher content of deoxyhemoglobin increases the rate of depolarization of hydrogen nuclei creating the NMR signal thus decreasing the intensity of the T2* image.The bottom line is that the intensity of images increases with the increase of brain activation. The problem remain that this increase is small (usually less than 2%) and easily obscured by noise and different artifacts.

2.1. Typical NMR imaging parameters

We use two different types of NMR images. The conventional images are obtain using our 1.5 Tesla GE 'Signa' scanner. Those images have good quality but it takes a few minutes to collect one, so they are used only to provide anatomical images of the scanned brain. The typical imaging parameters for the anatomical (T1 weighted) images are: FOV(Field of View 40 x 40 cm, slice thickness 7 mm, TE (echo time) - 11 ms, TR (repetition time) - 500ms, resolution 256x256 pixels. The functional images are obtained in EPI (Echo Planar Imaging) gradient echo sequence using the additional hardware from Advanced NMR. Those images can be obtained at the rate of several images per second. The usual imaging parameters are: FOV (Field Of View) 20*40 cm, slice thickness 7 mm, TE (echo time) - 45 ms, TR (repetition time) - 500ms to 3000ms, flip angle 60¡, resolution 128x64 pixels.

2.2. Main (known) sources of noise and artifacts

2.2.1. Thermal Noise

There is quite large thermal noise created in the electronic circuits of the MRI scanner. This noise can be quite well described as white uncorrellated noise. Its intensity-SNR (Signal to Noise Ratio) can be controlled only by increasing the pixel size thus loosing the spatial resolution.

2.2.2. Ghosts

Ghosts (strong artefacts looking like weaker copies of the actual images appearing next to it) may appear and disappear during one series of scans can create substantial spurious activation's. All series of images should be checked for such changing ghosts. Image series with ghosts should be discarded.

2.2.3. Cardiac and respiratory artifacts

The pulsation of the blood due to the heart rate and changes connected to respiratory cycle can change the blood flow and oxygenation. Those changes should have frequency higher than the changes in the activation pattern so that they will not create spurious activation's.

2.2.4. Draining Veins

The question whether the changes of the blood activation's come from the brain tissue itself or from the vein draining the active region is still open. Currently we do not use any additional steps to discriminate between those two contributions. There are two techniques that could be used. One calls for different imaging sequences (using spin echo or diffusion weighted imaging) that will be differentially sensitive to vessels of different size. Another approach is to use the fact that the blood in the draining veins is changing significantly (few seconds) later than the blood in the active tissue. So the activation's with longer delay can be interpreted as draining vessels. The most obvious precaution is to treat activation occurring close to visible large vessels with additional dose of skepticism.

2.2.5. Subject Motion

Subject motion is a serious source of artifacts. Even relatively small motion (of the range much smaller than a pixel size e.a. 1.6-3.2 mm) can create serious artifacts due to the partial volume effects (when the intensity of neighboring pixels varies by 10 % , which is quite common, then the shift of 0.1 pixel will create change of intensity of 1% - comparable to what we observe in the BOLD effect. Even small movements may also change the intensity of imeges due to the changes in spin magnetization during scan cycle. Subjects should be strictly instructed and constantly reminded not to move, straps restraing their head should be used. The task design should also minimize the possibility of task related movements. Since most of movements appear between consecutive series of images, only images taken within the same series should be compared directly in the data processing. We (mostly Chris Gatenby) are working on adapting the motion correction module to this program. This will be an additional program that will have to be run on the raw data before any further analysis.

2.2.6. Scanner Drift

There seem to be significant artifacts coming from the drift of the images in the readout direction. This drift is created most probably by the small instability of scanner gradients. It can create additional spurious artifacts in the regions that have large gradients of intensity in this directions.

2.2.7. Susceptibility artifacts (in EPI)

The EPI images are very sensitive to the changes of the magnetic susceptibility. In effect the signal from regions close to sinuses and bottom of the brain may disappear.

2.3.ÝReview of methods of statistical analysis

2.4. Advises for experimental design

Subject has to be repeatedly instructed that it is crucial for the fMRI study that he does not move at all.á Subjects responses has to be given in a way that minimizes motion - no speaking, only small movements etc.
Since the largest movements occur between imaging series the study should be organized so that tasks and baselines to be directly compared should be performed in the same imaging series.
To avoid drift artifacts (slow change of the overall image intensity over long time) each imaging series should start and end on the same task (ABABA, as oppose to ABAB).á First few images (usually 3 for TR = 1500) should be discarded because they are obtained before the machine reaches equilibrium. For each change of tasks the chemodynamical changes creating the change of intensity may take 3 to 8 seconds, so that first images obtained after each change should not be used.
Each imaging series should be repeated at least twice with the same task pattern to produce enough images to average.

Ý Back to the Table of Contents

3. Structure of this program

The core of this programming package is the function read_fmri that loads data from series of images and let user create and view different statistical images from this data. This can be used for a thorough interactive analysis of functional MRI data. fmri_batch is a version of this program that can be run as a batch job saving all prescribed statistical maps. To run read_fmri one have to create the setup file describing the structure of experimental data. This is done with help of functions try_setup, check_pos, hist_plot, fix_pos that help to asses the quality of the data file and define the appropriate strategy. After initial analysis the statistical parametric maps can be saved. For this purpose read_fmri can be run as a batch job, using its version called fmri_batch.

ViewRes provides access to all the statistical parametric maps created by read_fmri or fmri_batch and stored in our data base. It can display them using different threshold and different clustering filters. The interactive interface of ViewRes is very similar to that of read_fmri, but the former can acces only previously calculated statistical maps - all of them, while the latter is loading the whole raw data for a single slice of singe study - this is time consuming, but give acces to actual time course of intensities and can help to understand data in its complexity.

ÝResults of similar studies may be transformed into an uniform talairach space to be combined into composite statistical maps, or may produce sets of measures of activation's for predefined Regions of Interest (ROI).

Ý Back to the Table of Contents

4. Typical steps in the analysis of a functional project

4.1. Performed for each subject (imaging session)

Create a setup file (The easiest way is to copy setup file from a similar study and modify it).
Using try_setup and fix_pos adjust and test setup file
Analyze subject motion, drift and ghost arifacts, exclude bad images.
Run a batch job to create appropriate statistical maps.
Load data interactively (using read_fmri) and analyze time courses and statistical properties of interesting activated regions.
Print interesting activation maps.
Define Talairach coordinates.

4.2. Performed for the whole study

Create combined statistical maps for the whole study.
Define Region of Interest and obtain activation measures for each task, subject, activation task . This can further analysied using statistical software like SuperANOVA it should contains some error messages.

Back to the Table of Contents

5. How to get started on SUN or SGI workstation

5.1. login as the appropriate user (ask if you don't know the name

The program for the functional brain imaging fmri is available from every workstation as a user, using appropriate password.

5.2. start MATLAB

When in the proper directory start MATLAB by typing matlab in the command window. MATLAB is a program package (sort of programming language) which was used to write fmri. If you are running fMRI from other directory (your own or from the directory containing your data files) you have to copy the startup.m file to your directory. This will give matlab access to directories containing fMRI files.

Ý Back to the Table of Contents

6. Setup File - description of the study

The first thing you will need is to create the setup file describing all the parameters necessary in the data processing. This file should be located in the setups_fmri subdirectory of your working directory, its name should start from setup_fmri. Our usual convention is to use names like setup_fmri.7657od where 7657 is the study number and od characterizes the study type. The easiest way to start is to copy the setup_fmri.example file into your file here (named for example 9999X) and modifying it later. To copy this file type:

cp setups_fmri/setup_fmri.example setups_fmri/setup_fmri.9999XOther examples can be found at the end of this page

6.1. The File Structure of Data

The files created by the MRI are stored initially in the /anmr2/XYZ or /home/anmr/XYZ directory in the subdirectory called by the study number e.g.. 09999 (XYZ stands for initials of person who was performing scan). This directory has subdirectories (called series) named e.g.. 002, 005. Each series contains data from several experiments done with the same machine parameters setups. Within each series each consecutive run of the experiment consists in the sudirectory called image0, image1, image2, etc. Finally each of imageX directories contains actual image files . The data from consecutive slices are arranged in time order (If you have 3 slices and 20 images per slice taken during one minute scan, then images 0..19 belongs to the first slice, 20..39 to the second and 40..59 to the third slice. Images of all the slices are taken simultaneously, so that image 0, 20, 40 represents different slices of the image taken first, followed by 1, 21, 41).

6.2. Task Definition

To perform statistics you have to define which images belong to what activation. All images taken during the same activation belong to the same task. In the setup file you have to define how many different task were present and assign images to different tasks. Each image should be assigned to at most one task. Later tasks can be combined together. One task may contain data taken in consecutive runs (image0, image1,...). If you have changed the activation during scanning then you have to find out which images belongs to each activation. Usually you should skip about 5 initial images in each run (this is because the signal intensity is going down before the spin will reach equilibrium). If tasks are changed within one imaging sequence it is good to skip about 3 images after the task switch to account for the delay in the activation. (The chemodynamical equilibrium is believed to be reached after 2-8 seconds after onset of the task. The images that does not belong to any task are not loaded. If you want perform any data analysis on those images e.g. see the full time course of the activation, you may define additional task containing images not included in activation tasks (see the setup_fmri.example) those images may be then omitted while calculating statistics , but they images will be loaded and so represented in time course plots.

6.3. Line by line setup description of the setup_fmri.example

The # signs in the setup file are used as markers for the program that is reading it, so you cannot add or remove those signs (except adding or removing whole blocks or lines as explained later).

Study Name - any name you wish to give to your study
Study Number - if your study is on the optical disk give just its number (05835 in this case) you may also give the whole path of access beginning from /. Here the data is located in the /images directory.
Series Number - the number of series containing your functional data.
anmr directory name - usually initials of the person performing scan, the name of the directory in which the data where put while scanning.
preprocessing command - this command will be executed before reading data, usually left empty, this option can help you change some parameters during the run of the program.
Number of slices
Number of different tasks - number of different activation tasks performed. Those tasks (not necessary real tasks) defines the groups of images which will be treated together in performing statistics. Each image should belong to at most one task defined directly in the setup file. New tasks made from combinations of 'basic' tasks can be added using add_task postprocessing function. This number includes only tasks defined directly in the setup file - not tasks created using add_task.
The number of groups of lines following must be equal to the number of different tasks/sets defined here.
List of images for set this is list of subdirectories (called image0, image1, ...) containing the data belonging to this task. This must be followed by ranges lines, one for each image.
ranges - the image numbers from the appropriate image which belongs to this task. In this case we will have images ranging from 6 to 19 (inclusive) and from 46 to 59 belonging to each OFF E45 task and images 25 to 40 belonging to ON E45 task. There must be always even number of number between the # signs in each ranges line.
Series/File for anatomical - you may use anatomical images obtained using echoplaner (EPI), SIGNA or no anatomical images at all. If you leave this space blank the fmri will use the mean of the functional images as an anatomical image.
anatomical header size - 7168 for SIGNA and 40 for EPI image
anatomical image size - pixel number in X and Y direction (for SIGNA 256)
functional header size - header size for functional data (40 for EPI, 7168 for SIGNA)
functional image size (X/Y) - pixel size of the functional image
coordinates of interesting Field of View (FOV). Usually the image is much bigger than the size of the brain. To speed the processing and save the memory you should define the smaller rectangle containing only the important area. Here you give the coordinates of the upper left and lower right FOV corner. If you have a large study (100 images or more) it is really crucial to have a small FOV so that you can fit all the images in the computer memory. You should use try_setup to fine tune the size of FOV.
shift_x, shift_y - if you use the SIGNA anatomical image and EPI functional images they are sometimes shifted one with respect to another. After running the read_fmri program or from try_setup use fix_pos button to overlay those images. After finding the shifts values which seems to be best enter them to your setup file.
FOV factor (X and Y) - these define the ratio between the functional and anatomical machine defined field of view. They should be set to 1 if functional and anatomical are obtained in the same technique. In the typical EPI - SIGNA setup with the 20x40 EPI FOV and 40 x 40 SIGNA FOV, the FOV factor in the phase direction should be put to 2. This is usually the direction with the lower image size (Y in our example).
marker position - you may define the position of external water markers to load the data of their intensities (even if they are out of your chosen FOV). If you will not use it is better to leave this space free (this will save memory by not loading additional data).This is rarely used.
postprocessing commands - you may enter here the predefined postprocessing commands or any MATLAB commands which will be executed just after loading the data file. The most commonly used is the clean_data and init_med_filter command. For a description of available postprocessing commands go to Chapter 14.

6.4. check_pos : test for subject movements

To check your study for the subject movements you have to press the check_pos button. After a (looong) while you will see the new window with 5 plots. The upper two plots present the movements of the X and Y coordinates of the center of mass of the image intensity versus the image number (this represents time). The vertical axis is scaled in pixels. Usually if its variation is bigger than 1 pixel it means that there was quite a lot of motion. The third plot represents the number of pixel within this slice of the brain. It is designed to look for the movements in the z-direction. The forth plot presents the difference between consecutive images and the first image (solid line) and the last image (dotted line).Any big discontinuities in those plots usually indicates movements. Sometimes you can correct for it by eliminating the shifted parts of the data. The final plot presents the average intensity of signal from the brain used to normalize the signal intensity.

6.5. Ghost detection

The function check_pos calculate two factors that may help to detect ghosts: ghost and peak factor. Peak factor is sensitive to sudden differences ocuring in single images. Those images should be exluded from the analysis by adjustment in the setup file. The ghost factor is sensitive to the changes in the distribution of the intensity of pixels within the image. If either of them is larger than 10 it warns about the possibility of blocks of images with ghosts - ghost factor or about possibility of single images that differs suspiciously from neighboring ones. Closer study of histogram plot given by the function hist_plot can help to decide about the nature of the problem and possible exclude bad images from the setup file. Ghosts can be also observed in the window presenting currently loaded images.

6.6. SIGNA - EPI coregistration

The anatomical images from SIGNA sometimes are sometime shifted with respect to Echo-planar functional images. To use them as anatomical background for activation's they have to be coregistered. This is done by shifting images in plane by parameters shift_x and shift_y defined in the setup file. To find shift_x and shift_y one has to use the fix_pos function. It can be executed while running try_setup, or from the main program. The fix_pos window presents gray anatomical images with the edges of the first functional image overlaid in red. By changing shift_x and shift_y parameters you have to try to coregister as much anatomical structures of both images as possible. You may change the amount of edges displayed by changing the level parameter. The anatomical and functional images can be also seen next to each other or one above the other. The best shift_x and shift_y values should be entered to the setup file.

6.7. try_setup : previewing data

try_setup is a function for trying the setup file without loading all the data sets. It presents the anatomical images for all the slices and overlays edges of functional images. It should be always run after creating the setup_file to ensure that it works.

6.8 example of setup file.

never add or remove any '#' sign from this file #
Study NameÝÝÝ # dyslexia project6 7-24-96 Lee Chupka #
Study NumberÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ # /iguana1/project6/10050 #
series numberÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ # 004 #
anmr directory nameÝÝÝÝÝÝÝÝÝÝÝÝ # RTC #
preprocessing commandÝÝÝÝÝÝÝÝÝÝ # #
Number of slicesÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ # 8
Number of images per sliceÝÝÝÝÝ # 30
Number of different setsÝÝÝÝÝÝÝ # 5

List of images for set # 0 5 10 15 #

OFF set name # line #

ranges # 4 30 # ranges # 4 30 #

ÝList of images for set # 1 6 11 16 #
ON set name # case #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #

ÝList of images for set # 2 7 12 #
OFF set name # letter rhyme #
ranges # 4 30 #
ranges # 4 30 # ranges # 4 30 #

ÝList of images for set # 4 9 14 19 #
ON set name # category #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #

ÝList of images for set # 3 8 13 18 #
OFF set name # non-word rhyme #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #
ranges # 4 30 #

Ý# end_of_range_lists #

ÝSeries/File Name for anatomical # 002 #
anatomical header size # 7168
anatomical image size -Y- # 256
anatomical image size -X- # 256

Ýfunctional header size # 40
functional image size -Y- # 64
functional image size -X- # 128

ÝCoordinates of interesting FOV
upper left cornerÝ
y = # 2
x = # 42
lower right corner
y = # 63
x = # 89
shift_x # -1
shift_y # 1
FOV factor - X # 1
FOV factor - Y # 2

Ýmarker position # #

Ýpostprocessing comands :
# clean_data #
# init_med_filt #
# end_of_file #

Ý Back to the Table of Contents

7. Running read_fmri

Starting program After creating the setup_fmri.9999X file in the setups_fmri directory you have to be back in your main directory /home/oursun/fmri. At first you have to start MATLAB by typing matlabÝ. The square matlab window with a funny graph should appear for a while and disappear, your prompt character will be now >>. To start program you type read_fmri you will see the list of available setup files and you have to chose one and type its name (only the part after the dot e.g.. 9999X). Later you have to define the slice number. You can process only one slice at time. The program will start loading consecutive image files printing its names. It will also display the small images titled by the image numbers (1_34 is the image number 34 from the subdirectory image1). You should watch those images to check whether there are no huge artifacts suddenly appearing in certain images. If this happens you should modify your setup file so that you don't process those spoiled images. After successful loading of the data (it may take couple of minutes) the newimageEPI window with plenty of buttons and with your image should appear.

8.1. Choosing tasks

The long ON and OFF buttons in the lower part of the main imageEPI window are defining the sets(tasks) which will be used to the current statistics. After pressing them you can chose between task as they were predefined in the setup file. You can also use (All) images, or (for ON button) (All but OFF) to define all images that are not included in the current OFF set.

8.2. Choosing statistics

By pressing the statistics button you can chose from a wide variety of statistics to be performed on the currents task. The statistics will be performed for each pixel separately on the sets of intensities found for this pixel in the images belonging to the currently chosen ON and OFF sets. This will create a map of intensities (called here t-map even if created with some other statistics). To actually perform the statistics after changing the statistic or ON, OFF sets, filters definitions press run_stat. Here are most commonly used statisctics from the list of statistics that you can choose from :

Ýdirect t-statistics % 1
Ýspilt-2 t-statistics 2
Ýcombined t-statistics % 5
Ýpercent difference % 6
combined (using mean) skew statistics % 26
combined skew statistics % 27

check also the review of statistics in chapter 2.3

8.2.1.ÝSplit Statistics

The split statistics are statistics performed on parts of data separately and added logically, e.g. for split-2 Mann-Whitney the ON and OFF sets are divided into two equal parts, Mann-Whitney statistics is performed on both subsets, and two 't-maps' are created. Only activation that are stronger than the chosen t-cutoff on both maps are displayed as positive activation's. Those statistics seems to be batter if your data are composed from several blocks, performing statistics on each block separately may be better to avoid motion and signal drift artifacts. Split statistics divides On and OFF blocks into equal parts, so it can be used only if the blocks of images obtained in different imaging series are exactly the same length in other case (e.g. if some images were removed due to ghosting) the combined statistics have to be used.

8.2.2.ÝCombined Statistics

The combined statistics is an improved version of split statistic. It performs statistics between images that belongs to the same imaging series. For each imaging series it checks if there are both ON and OFF images in this series. If it is true, it calculates the statistical map from images found in this series. The maps calculated for all the series containing both ON and OFF images are later combined. If there are less than for such maps all have to be larger than the cutoff, if there are between 4 and 7 maps all but one has to be larger in general Int(n/4) of maps are allowed to be smaller than cutoff. If the study contains equal number of ON and OFF images in each imaging series then combined statistics is exactly equivalent to the split statistics with spliting factor equal to number of imaging series containing data.

8.2.3.ÝSkew Statistics

The skewed t-statistics is similar to t-statistics, but it takes into account the linear drift that may apperear in the data. This is happening quite often mostly due to the slow machine drift that shifts images over a fraction of pixel during the imaging series. As the regular t-statistics approximates data by a function f(t) = a*t(ON) + b*t(OFF), the skew statistics approximates it as: f(t) = a*t(ON) + b*t(OFF + c*t. The power of difference is calculated as difference a-b divided by the standard deviation of the data from this fitting function. This statistics is really helpful if the check_pos show that a small and consistent drift exist in the data.

8.2.4. ROIcorrelation

The ROI correlation require that certain Region of Interest (ROI) is specified (see ROIinterrogation). This statistics calculates the correlation coefficient between the pixel time course and the average time course of the chosen ROI.

8.2.5. AutoCorrelation

- widths of autocorrelation calculates the smallest lag at which the autocorrelation falls below 1/3 of its value at the zero lag. Large width of autocorrelation means long temporal correlation present in the data. task autocorrelation the autocorrelation on the lag defined by the task switching pattern normalized by the zero lag autocorrelation (between 1 and -1).

8.2.6. Fourier Statistics

This statistics calculates the power of specified Fourier components. In viewing results one can choose the phase window and frequency of stimulus.

8.3. Cutoffs and filters

The controls for the cutoff value and filters are located in the bottom left corner of the imageEPI window. Each pixel which have the current statistics value bigger than the cutoff will be displayed in hot (red or yellow) color, pixels below negative cutoff are displayed in cold (blue or violet) color. In the upper corners the color maps and their limiting values are displayed. The numbers below describe the number or positively and negatively activated pixels. By changing the cutoff value you may change the significance level, changing the number of activated pixels. To actually perform the statistics after changing the statistic, ON, OFF sets, cluster filters or t_cutoff definitions press run_stat.

8.4.ÝCluster and neigborhood filters

The spatial size of activated are is usually larger than one pixel. In other words the spatial correlation of activation map is much larger than spatial correlation of noise. This can be used to increase power of our statistical analysis by accepting as real activations only pixels that belongs to the larger activated area. This can be done by applying certain cluster or neighborhood filter to the thresholded activation map. Those filters can be aplied by choosing the parameter cluster_number different than zero. This can be done while viewing activation maps using read_fmr, or ViewRes.

8.4.1.ÝNeighborhood filter

The positive value of cluster_number means that the neighborhood filter will be applied. For each activated pixel this filter counts the number af activated neighboors (wall neighboor counts as 2, corner neighboor as 1). The pixel is presented as activated only if the sum of neighboors is larger or equal cluster_number. This filter can still leave isolated pixels (when pixel had many neighbours, but those pixels did not have enough). Basically it strips the outside pixels from each activated cluster. After applying this filer even single surviving pixel represents larger activation.

8.4.2.ÝCluster filter

If cluster_number is negative it means than another cluster filter will be applied. For each pixel the size of cluster that it belongs to decides whether it is treated as activated. Only pixels belonging to clusters greater than -cluster_number are presented as activated. Unfortunately this filter cannot be applied if the threshold is too low e.a. if substantial part of the whole image is above threshold.The advantage of this filter is that it does not take out pixels that belongs to extremities of large clusters. It also makes clearer images because it does not leave small cluster or isolated pixels.Ý Back to the Table of Contents

8. Additional data analysis

9.1. ROI interrogation

The time course of the signal intensity for smaller region of interests ROI or even single pixels may be displayed using the buttons located above the image. The ROIbox button will define the rectangular region with two corners defined by mouse. The ROI act+ (ROI-) will take only the positively (negatively) activated pixels from the marked rectangle. The upper plot represents the signal intensity averaged within defined ROI as a function of image number (time). The lower plots overlays time-course curves for all the pixels within chosen ROI (shifted by their average intensity). Labels give the average signal intensity for ON, OFF images and the average values of the current statistic. The bottom histogram represents the distribution of current statistic between ROI pixels.

9.2. ROIedit : more precise but tedious ROI definition in MATLAB

By clicking on ROIedit you can open the ROI editor that will help you define the more complicated ROI, by adding and deleting single pixels or boxes. You can also save and load ROI's. The ROI will be stored in the directory for saved data defined in your setup file (default is /home/mrpart/fmri) in the subdirectory ROI-files and further in the subdirectory named after your setup file name. If you select the add_box or delete_box button in the ROI editor you will have to move your cursor to the image (it will change into a small cross)and click in two opposite corners of the selected box. If you select add_point or delete_point you will be able to click on individual points, the changes will be introduced after you will hitÝto execute the addition (or deletion).If you want to define ROI's for future use in a study in which you have already saved some data (anatomical images and statistics maps) using batch jobs or save_data, you don't have to load the whole functional data. You can view anatomical images and use them to define ROI using the view_anat function. Just type it in your command window (with MATLAB running).

9.3. ROI definition using Adobe Photoshop

Defining the Regions of Interest by drawing them on anatomical images is very time consuming. This can be speeded up by using drawing abilities of Adobe Photoshop. The anatomical images have to be saved in TIFF format using MakeTiffFolder those images created in the /oursun1/images directory have toe transferred using fetch to Macintosh that can run Adobe Photoshop. Those anatomical images uses colors with numbers from 11 to 210. The other colors can be used to paint regions of interests. Anatomical images with ROI painted over should be transferred to the main data base and can be later used by the ROIoutput function to generate measures of activation's for those regions.

9.4. save_data : Saving t-maps

You can save your current t-map and anatomical image by hitting save_data. The data will be saved in the same format as in the batch job, they can be wieved using tools for viewing batch job results ViewRes.

9.5. Markers

If water marker were used during the experiments you can utilize them to obtain the time course of signal from those markers as well as motion of their center of mass. To use them you have do define its position in the setup file. You have to give two coordinates of the center of marker and you may also give the third parameter s defining size of the marker window to be 2*s+1 (default is s = 3). If your setup file contains marker coordinates you will have marker_plot button in the imageEPI window. Hit it to obtain the time course of the marker center of mass, its mean intensity, as well as an image of marker (to check if its coordinates are correct). To find coordinates of marker you can do following. After having loaded your data using the setup file without marker coordinates type figure in the matlab window. Then type image(im_data1') this will create image of one of the full functional images. Next type ginput(1), move your cursor to the center of visible marker and hitÝ. The marker absolute coordinates that will be printed should be now copied into the setup file.

9.6. Viewing raw images

The raw images buton on the main interactive window opens window for displaying raw images in the currently loaded image series. The images can be called by number, played as a movie or displayed in an array of 25 equally spaced (in time) images. The raw intensity images can be viewed as well as difference images (with median or neighbooring image subtracted.

Ý Back to the Table of Contents

9. Correlation analysis

Statistics 11 - 16 are design for analysis based on the correlation coeficient between the pixel time course and activation time course.

ÝThere are three option to be used as a reference time_course:

The time_course of particular region of interest.
Box-car function defined to be 1 for ON images and -1 for OFF images
user defined time_course predefined in the global variable time_course

For each of those options ther are also two possibilities differing with the treatment of individual series:

whole time course of images belonging to ON and OFF task is taken into account - this may be dominated by the effect of changes in intensity between series.
correlation coefficient calculated for each series separately and mean of those coefficients used in the final output. In this case whole series containing images from ON and OFF task are taken into account.

Following functions perform correlation statistics

11. ROIcorel
Correlation to the timecourse of predefined Region of Interest (ROI)

12. ROIcorel_comb
Correlation to the timecourse of predefined Region of Interest (ROI) calculated for each series (containing images from ON or OFF task) separately and combined.

13. block_corel
Correlation to the block (box-car) function.

14. block_corel_comb
Correlation to the block (box-car) function calculated for each series (containing images from ON or OFF task) separately and combined.

15. ext_corel
Correlation to the user defined (external) time course defined as a global variable time_course.time_course should be an vector of length equal to the number of images.

Ý16. ext_corel
same as ext_corel, but calculated for each series separately and combined.

Ý Back to the Table of Contents

10. Fourier Analysis

Back to the Table of Contents

11. Running a Batch Job

You can run your statistics overnight by running a batch job. To do this you have to create your batch file. You can do it by copying the batch_example into your file. The parameters to be modified include:

Ýbatch_example.out (in the second line) - the file that will contains all theÝName : name of the setup file
Slice : list of slices on which the statistics will be calculated,
Stat : list of statistics to be calculated the numbers refer to the order they appear in the main imageEPI window when you are choosing statistics manually for example (1-split4T, 2 - split2T, 5 - T-statistics, 6 - % difference).
OffOn : is the list of pairs of tasks on which the statistics are to be performed the numbers follows the order in which tasks are defined in the setup file.

To run a batch job named my_batch.csh on 11PM you have to be in the command window and type:

at -c 2300 my_batch.csh (on oursun)
at -f my_batch.csh 2300 (on all other)
If you are still running MATLAB (have >> as a prompt) you may execute system commands preceding them by the exclamation mark !

!at -f my_batch.csh 2300 (on mistral, from MATLAB)
The batch job will create the directory named after the setup file in your 'save data directory' - default name is /fmri_data. separate file will be created for anatomical image of each slice named e.g. anat5835a_1.mat. and for t-maps e.g.res5385a_12sl1st2 (this means statistics number 2 performed between tasks 1 and 2 on slice 1).

Example of the batch file batch_example

Back to the Table of Contents

12. Viewing Saved Results

12.1. Viewing multiple maps

To view the results of the batch job (or t-maps saved using the save_data command) you use the ViewRes command. The format is following:

ViewRes(setup_name, pair_list, t_cut, cluster_nb, slice_list, stat_nr,options)

for example:ViewRes('9999X','[12 34],2.3,5,[1 2 3],2,'g0')will create a window with images saved for the study 9999X, presenting the comparisons between the task 1 and 2 in the first row and 3 and 4 in the second row using statistic number 2 (this is split 2 T stat.). The cutoff parameter will be 2.3, cluster filter of 5 will be applied. Options : '0' and 'g' will be appled to the image to use the t_tab0 map (instead of default t_tab1) and to gausian filter the final map.

To print this image type print_image in the command window. This will transform color activation dots into white (for positive) and black (for negative activation's) and print this image. To create more customized you can make your copy of /fmri_data/fMRI-view/ViewResMacro.m file and change it .

12.2. Interactive viewing

By typing ViewRes without parameters one can invoke the interactive window (similat to the main fmri_read window to see any statistical map that is strored in the data base. If study have composite maps obtained using Tal Compose or TalCompose3D

Ý Back to the Table of Contents

13. Registration of images into the Talairach Atlas

To combine data taken from different subjects and to report them in terms of Talairach -Tarnoux coordinates the images may be transformed into the standard from. To make it possible the images have to be taken in a plane exactly perpendicular (coronal) or parallel (axial) to the AC-PC line connecting the anterior commisure with posterior commisure. After defining necessary talairach coordinates each image may be transformed into the uniform talairach space. This operation can be also performed for the whole stack of images 3-dimensional talairach transformation.

13.1 In plane talairach definintion.

First one have to run the function tal_define to define the direction on the midline and positions of AC, PC and the furthest extent of the brain. For the setup name 7777x it will create file named tal7777x in the /fmri_data/7777x directory. < H2>13.2 Definitions for 3-dimensional transformation.If the positions of slices with respect to the AC-PC line differs between subjects one have to add information from the sagital image to the file created by this function. This has to be done always when more then 3 slices are taken. To find those information you need to have a printout of the sagital image with lines marking the position of slices. If the sagital information are given the 3D version of programs using those coordinates have to be used later.

13.3 Format of talairach space.

After the Talairach coordinates are defined they can be used to transform anatomical images and statistical maps into the standard talairach form. For axial images it will be array of size (40+15+40)x(40+40) pixels, for coronal (40+20)x(40+40) pixels. Those images may help in finding the talairach coordinates of pixels, they may be also used to make combined images from several studies and to calculates statistics for regions of Interests defined uniformly in the Talairach coordinates system.

13.4 Studies taken on an angle to AC-PC system

For studies in which slices are taken in the different plane the talairach functions may be used for creating composite images or defining Regions of Interest by creating "talairach-like" coordinate system based on different reproducible anatimical landmarks.

Ý Back to the Table of Contents

14. Combining data from several studies

Quite often the changes created by the differences in tasks are so weak that for any single subjects we are not able to find activation's that are significantly stronger than random false activation's. In some other studies activation's vary a lot between subjects. These some of the reasons to use the data from several subjects performing the same tasks to increase the statistical power of our analysis. This procedure used always in PET studies. Early fMRI studies of visual and motor cortex showed enough signal to pinpoint activation's in single subjects. Yet when the focus shifted to image more subtle cognitive functions, certain techniques of combining data from larger groups of studies have to be used. Currently we use three methods of analyzing data between subjects.

14.1. ROI analysis

The most popular way is to define Regions of Interest (ROI), calculate the measure of activation for each ROI, each task and each subject and perform ANOVA or other statistical analysis on this reduced data set. There are two ways of defining ROIs. To give a precise anatomical localization they can be drawn on each anatomical image using tools of ROIedit or Adobe Photoshop for each subject separately. The other less time consuming way is to define them once in the Talairach coordinate space and apply Talairach transformation on the data to find actual position of ROI for each subject. Due to differences on individual anatomy this method can be used only with large ROI. Usually they are defined as unit blocks of Talairach coordinate system or at most half of this blocks (5x5) pixels in our resolution of talairach transformed images. To define set of ROIs in Talairach space use ROI3Ddef. To create file with measures of activity in the regions that can be later processed using Super ANOVA, EXCEL or other programs use DoRegions.

14.1.1 Measures of ROI activity

The following measures describing the activation in the ROi can be calculated:

n - number of active pixels (defined as pixels above certain threshold and satisfying certain cluster filter)
p - mean value of statistics for the whole region - mainly to use with percent difference.
s - sum of values of statistics 1 (usually percent difference) summed only for pixels recognized as active by statistics 2
r = s/n
q - number of active pixels normalized by the region size
Q - number of active pixels normalized by the region size and total number of activation's in the whole brain
x,y,z - coordinates of the center of intensity for activated pixels
R - dispersion of intensity for activated pixels (mean distance form the center of mass.

14.2. Combination of statistical maps

Data form several statistical maps may be combined together using the TalCompose or TalCompose3D functions. The latter have to be used if slices are not equivalent and sagital information are entered while defining the Talairach coordinates. The composite anatomical images can be made by use of TalComposeAnat or TalComposeAnat3D

14.3. Uniform t-statistics on data sets combined from many subjects

T-statistics may performed treating all the images from different subjects (transformed into uniform Talairach space) as one uniform set of images. To do this one have to calculate statistical maps for each subject using statistics 17. The data files created for this statistic will contain enough information to calculate the final t-map. After all the studies for different subjects have Talairach coordinates defined and statistical maps st17 are created the function TalTotalStat have to be executed to create the combined statistical map.

Ý Back to the Table of Contents

15. Postprocessing Commands

The following postprocessing commands can be executed inter actively after loading data (using read_fmri or fmri_batch). They can be also added to the setup_file (in the end of the file between # # signs), in this case they will be executed just after loading data by either read_fmri or fmri_batch.

add_task

You may create a new task by adding already defined tasks by adding this command to the setup file as a postprocessing command. The format is: # add_task([1 3 4],'combined task') #This will define a new task being a sum of earlier in this setup file defined tasks 1,2,4 and called 'combined task'. This task will be assigned first available number. This command may be applied also from the MATLAB command line (without # signs). This command should be always used when you want do define different tasks containing the same images.

init_med_filt(s)

This will perform the median filter on each image (replacing each pixel intensity by the median of the s x s neighborhood, default value is 3), For large data sets it is a time consuming option. The median filter is applied to the difference between each image and the mean image, to do this the mean image is first subtracted from each image, median filter is performed and the mean image is added back.

init_gaus_filt(s)

Performs a gaussian filter of width s on all the data images.

threshold(f)

the background is cleaned by zeroing each pixel with intensity below f of the maximum intensity (default value is 1/12). This helps to remove outside noise and to calculate center of mass of the image intensity. If you are interested in the area where signal was surpressed by susceptibility effect you should lower the threshold value. It is part of clean_data postprocessing command.

time_normalization

the image intensity is normalized in time so that the average intensity of the whole image is constant. Each image is divided by its average intensity and multiplied my the average intensity of the whole study. The lowest plot in the check_pos and marker_plot windows will present the time course of this normalization constant.

clean_data

this option is used most often it combines the threshold and time_normalization commands.

spatial_filtering

the possibility of spatial filtering to eliminate the low frequency spatial signal changes are being considered - still on the experimental level.

Back to the Table of Contents

16. Dictionary of useful commands

BackupData, Restore data Moves whole directory from the /fmri_data to different location, checks

batch_fmri

Runs the main program as a batch job for specified slices, statistics, and task comparisons

check_pos

creates series of plots for analysis of movement

DoRegions

Calculates chosen measures of activations (number of active pixels, percent difference of activity or their combinations) for whole group of setup files and create text files for analysis with other software (super ANOVA).

fix_pos

helps adjust shift_xb and shift_y to corregister anatomical and functional images.

FourierCompose, FourierDiffer

tool for displaing the Fourier analysed data combined from several imaging series. Displays pixels that has common phase or that differs by defined phase lag (e.g. when order of activation tasks is switched).

hist_plot

plots of distribution of pixel intensities in time - tool for analysing the ghosts artifacts.

imageEPI

runs the interactive window for data analysis (run automatically by read_fmri).

MakeTalFolder

MakeTiffFolder

Create a folder with anatomical images in the TIFF format in the /home/hundun_mac directory (/home/mac as seen from hundun). This helps to move images into Adobe Photoshop on the Mac to draw ROIs by hand.

print_image(orientation)

changes the color map from into bw_col_map containing white for all the posiive activations and black for negative ones and prints the current window on black and white laser printer. orientation may be defined as 'tall', or 'landscape'.

RandomDataSet

Creates random data set - for testing new functions.

read_fmri

runs the main program to load the images and analyse them in the interactive mode.

ROIoutput

ROIprint

outputs the intensities of pixels in the chosen ROI.

ROI3Ddef

tool for defining ROIs composed from blocks in the talairach space.

talairach

transforms image (anatomcal or statistical map) into the uniform talairach space - one slice only

talairach3D

transforms 3D set of images of several slices (anatomcal or statistical map) into the uniform talairach space

TalCompose,ÝTalComposeAnat

Creates a composite activation map by transforming images into the Talairach space and averaging them (using mean or median). Assumes that slices in each study are equivalent.

TalCompose3D,ÝTalComposeAnat3D

Creates a composite activation map by transforming images into the Talairach space and averaging them (using mean or median).

TalDefine

tool for definig the Talairach coordinatesin anatomical images.

TalTest

Checks the accuracy of the Talairach definition - presets the anatomical images transformed into the Talairach space. (in slice only)

TalTest3D

Checks the accuracy of the Talairach definition - presets the anatomical images transformed into the Talairach space using (3D transformation)

TalTotalStat

try_setup

tool for testing the setup_file. displays anatomical images and overlays functional images in them.

ViewRes

without parameters, runs interactive tool for displaing saved results.
Ýwith parameters displays results for specified setup_files, task comparisons,slice, cutoff, cluster filter

ViewResTal

displays saved results for a study transformed into th talairach space.Ý Back to the Table of Contents

17. Using optical disks

17.1. Small (128 MB) optical drives

to mount : mopt
to dismount : dopt
to format : format_optical.
the mounted disk directory is : /usr/optical

17.2. Large (1.3 GB) optical disks

to mount : opcom mbigopt
to dismount : opcom ubigopt
the mounted disk directory is : /usr/bigopt
Large optical disks are installed on Oursun, Mistral, Boreas and Iguana.

attention: to dismount disk you must be sure it is not busy. It means nobody (in no window) is in the directory of this optical disk

Ý Back to the Table of Contents

18. Some technical details

The files necessary to run fmri are located in the following directories in /fmri_data :

fMRI-tools - main toolbox
fMRI-post - postprocessing commands
fMRI-view - commands to view saved statistical maps

All the global variable necessary to run fmri are defined and described shortly in the file globals.m in fMRI-tools.Ý

19. Troubleshooting ( before you ask Pawel)

19.1. try_setup creates strange, doubled or distorted images

check the image size , Field of Viev (FOV), and sscaling factors, try switching x and y directions.

19.2. errors appearing while reading setup file

do all the # signs for strings have their closing partnersis the number of tasks validwithin each task there should be as many lines as many image directories contain this task.

19.3. out of memory error while loading data

Quite often program breaks at this moment with the out of memory error. This means that you are trying to load to many images into the memory. To help it you should:

check if somebody else is running huge jobs on your computer,
decrease your Region of Interest cutting it closer around the brain,
eliminate some activation's from your setup file (you may create different setup files to perform different comparisons),

19.4. your batch job didn't produce any files

check the batch.out file created by this job. If it doesn't exist it is most probable that the job never run. If it exist it should contains some error messages.

Ý Back to the Table of Contents

20. Short Description of Statistics

1) - regular t-statistics
2) - split 2 t-statistics
3) - split 3 t-statistics
4) - split 4 t-statistics
5) - combined t-statistics (combined between series)
6) - percent difference
7) - combined percent difference
8) - Mann-Whitney test
9) - combined Mann-Whitney test
10) - paired t-statistics
11) - ROI correlation
12) - ROI correlation combined
13) - OnOff Block correlation
14) - OnOff Block correlation combined
15) - correlation to externally defined function
16) - correlation to externally defined function combined
17) - data for performing t-stat for several studies
18) - mean intensity of ON images
19) - standart deviation of intensity of ON images
20) - mean difference with skew correction
21) - Fourier power spectrum
22) - percent difference with skew correction
23) - t-statistics with skew correction
24) - split 2 t-statistics with skew correction
26) - skew corrected t-stat (combined using mean)
27) - skew corrected t-stat (combined using minimum)
28) - slope used for the skew correction
30) - neighboors correlation
31) - gradient in X direction
32) - gradient in Y direction
33) - Fourier phase
34) - gamma1 - skewness
35) - gamma1 - kurtosity
45) - AutoCorrelation
46) - Talairach Block time courses
52) - randomized (bootstrap) statistics

ROI time course

Statistic 46 calculates the time course for chosen tasks averaged for the basic Talairach blocks ROItal_tab_off, ROItal_tab_on.
For tasks that contain series of blocks of exactly equal lengths the time course is averaged to one epoch.
The results can be viewed using functions ROItc_disp(ROItal_tab_on), or ROItc_disp2(ROItal_tab_on,ROItal_tab_off). Those variables may be calculated in batch job and later averaged between several studies.

ÝBack to the Table of Contents