Evaluation and implementation of neural brain
activity detection methods for fMRI
Sabina Breitenmoser
2005-02-21
´
Ecole Polytechnique F´ed´erale de Lausanne Swiss Federal Institute of Technology
School of Engineering - Signal Processing Institute (ITS)
EPF Lausanne , Prof. Jean-Philippe Thiran
Institute of Technology - Biomedical Engineering (IMT)
LiTH Link¨oping , Prof. Hans Knutsson
Link¨oping University, Universitetssjukhuset SE - 581 85 Link¨oping
Evaluation and implementation of
neural brain activity detection
methods for fMRI
Diploma Thesis
Winter Term 04/05
Professor EPFL:
Prof. Jean-Philippe THIRAN
Professor LiTH:
Prof. Hans KNUTSSON
Assistant:
Joakim RYDELL
Sabina BREITENMOSER
Section de g´
enie ´
electrique et ´
electronique
Avdelning, Institution
Division, Department
Linköpings tekniska högskola
Institutionen för medicinsk teknik
Datum Date 21 februari 2005 Språk Language Rapporttyp Report category ISBN Engelska/English
Examensarbete ISRN LiTH-IMT/ERASMUS - R - - 05/030 - - SE
Serietitel och serienummer
Title of series, numbering
ISSN
URL för elektronisk version
http://www.ep.liu.se/exjobb/imt/erasmus/2005/030/
Titel
Evalutation and implementation of neural brain activity detection methods for fMRI Författare
Author
Sabina Breitenmoser
Sammanfattning
Abstract
Functional Magnetic Resonance Imaging (fMRI) is a neuroimaging technique used to study brain functionality to enhance our understanding of the brain. This technique is based on MRI, a painless, non-invasive image acquisition method without harmful radiation. Small local blood oxygenation changes which are reflected as small intensity changes in the MR images are utilized to locate the active brain areas. Radio frequency pulses and a strong static magnetic field are used to measure the correlation between the physical changes in the brain and the mental functioning during the performance of cognitive tasks.
This master thesis presents approaches for the analysis of fMRI data. The constrained Canonical Correlation Analysis (CCA) which is able to exploit the spatio-temporal nature of an active area is presented and tested on real human fMRI data. The actual distribution of active brain voxels is not known in the case of real human data. To evaluate the performance of the diagnostic algorithms applied to real human data, a modified Receiver Operating Characteristics (modified ROC) which deals with this lack of knowledge is presented. The tests on real human data reveal the better detection efficiency with the constrained CCA algorithm.
A second aim of this thesis was to implement the promising technique of constrained CCA into the software environment SPM. To implement the constrained CCA algorithms into the fMRI part of SPM2, a toolbox containing Matlab functions has been programmed for the further use by neurological scientists. The new SPM functionalities to exploit the spatial extent of the active regions with CCA are presented and tested.
Nyckelord
Keyword
Contents
1 Introduction 11
2 Motivation 13
2.1 Interaction between engineering and medicine . . . 13
2.2 Medical signal processing and neurology . . . 13
3 Background 15 3.1 Magnetic Resonance Imaging (MRI) . . . 15
3.1.1 Introduction . . . 15
3.1.2 Basic principles . . . 15
3.1.3 Imaging parameters . . . 17
3.1.4 Field of application . . . 18
3.2 Functional Magnetic Resonance Imaging (fMRI) . . . 19
3.2.1 Introduction . . . 19
3.2.2 The Blood Oxygen Level Dependent (BOLD) response . . . 19
3.2.3 Brain activity detection . . . 20
3.3 Statistical Parametric Mapping (SPM) . . . 20
3.3.1 Introduction . . . 20
3.3.2 SPM approach . . . 20
3.3.3 SPM image format . . . 21
4 Theory 23 4.1 GLM method . . . 23
4.2 Canonical Correlation Analysis (CCA) . . . 24
4.2.1 The canonical correlation . . . 24
4.2.2 Constrained CCA . . . 26
4.3 Principal Component Analysis (PCA) . . . 26
4.4 Impulse response . . . 27
4.5 Spatial filter functions . . . 27
4.5.1 Gaussian filters . . . 27
4.5.2 Steerable filters . . . 28
4.6 Statistics . . . 30
4.6.1 ROC (synthetic datasets) . . . 30
5 Methods 33
5.1 Evaluation of neural activity detection methods . . . 34
5.1.1 Brain data acquisition . . . 34
5.1.2 Brain data handling . . . 34
5.1.3 Implementation of the CCA method . . . 35
5.1.4 Generation of a hemodynamic response model . . . 35
5.1.5 Spatial modelling . . . 38
5.1.6 Generation of a correlation map . . . 41
5.1.7 Noise filtering . . . 42 5.1.8 Statistics . . . 42 5.1.9 Matlab GUI . . . 43 5.2 SPM implementation . . . 44 5.2.1 Structure of SPM software . . . 44 5.2.2 Testing SPM2 code . . . 45
5.2.3 Integration of CCA methods . . . 45
6 Results 47 6.1 Some important values . . . 47
6.2 Evaluation on synthetic data . . . 47
6.2.1 Related work . . . 47
6.3 Evaluation on real human data . . . 48
6.3.1 Steerable spatial filters . . . 49
6.3.2 Constrained CCA . . . 49
6.4 Implementation in SPM . . . 50
6.4.1 New SPM code for CCA . . . 50
6.4.2 New SPM architecture . . . 52
6.4.3 New functionalities . . . 52
6.4.4 Evaluation in SPM: GLM versus constrained CCA . . . 52
7 Discussion 57 8 Conclusion 59 A Guideline to use SPM for fMRI 65 A.1 Getting started . . . 65
A.2 The SPM menu window . . . 65
A.3 A canonical correlation analysis with SPM . . . 67
Abbreviations
BOLD Blood Oxygen Level Dependent bpm beats per minute
CCA Canonical Correlation Analysis EPI Echo Planar Imaging
fMRI functional Magnetic Resonance Imaging FWHM Full Width at Half Maximum
GLM General Linear Model GUI Graphical User Interface
hrf hemodynamic response function ICA Independent Component Analysis MIP Maximum Intensity Projection MRI Magnetic Resonance Imaging NMR Nuclear Magnetic Resonance PCA Principal Component Analysis ROC Receiver Operating Characteristic SNR Signal to Noise Ratio
SPM Statistical Parameter Map TR Time to Repeat
Chapter 1
Introduction
The brain is the most fascinating and least understood organ in the human body. For centuries, scientists and philosophers have tried to find the relationship between behav-iour, emotion, memory, thought, consciousness and the physical body. In recent years, techniques for non-invasive monitoring of the working brain have experienced a strong development. Neuroimaging tools and methods have been designed to study brain func-tionality to enhance our understanding of the brain.
Functional Magnetic Resonance Imaging (fMRI) is one neuroimaging technique with the capacity to map neural activity with high spatial resolution. The technique is based on MRI, a painless, non-invasive image acquisition method without harmful radiation. In fMRI radio waves and a strong magnetic field are used to measure the correlation be-tween physical changes in the brain and the mental functioning during the performance of cognitive tasks. The physical changes are small intensity changes in MR images due to local blood oxygenation changes.
This project presents some approaches for the analysis of fMRI data. The Canonical Correlation Analysis (CCA) is evaluated and studied. CCA is able to fully exploit the spatio-temporal nature of fMRI data for detecting active brain areas. In a second step, an implementation of this promising technique into the software environment called SPM (Statistical Parameter Mapping) is performed. SPM is designed to analyse fMRI data as well as PET/SPECT images. The spatio-temporal CCA-based fMRI techniques are not previously supported by this software.
Chapter 2
Motivation
2.1
Interaction between engineering and medicine
Cross-disciplinary activities in the field of medicine are of a growing interest. Biomedical engineering is a discipline that joins the knowledge in engineering sciences with clinical practice. New devices, algorithms, processes and systems are developed to improve the medical practice and health care. A researcher in the field of biomedical engineering not only needs to be familiar with the relevant applications of engineering in medicine but also with the basic life sciences.
This interaction between the traditional engineering field and modern medicine is for me the motivation to tackle this project. It is interesting to see how new techniques in the field of engineering can improve medical diagnosis and health care.
2.2
Medical signal processing and neurology
Medical care incorporates diagnostic, monitoring and therapeutic issues. Typically rele-vant patients characteristics are measured, interpreted and an appropriate decision about the therapeutic action is taken. The medical signal processing focuses on the extraction of the useful information from medical images and signals.
In this project we will work with medical MRI images of the brain to extract neural brain activity. The interpretation of this information and the decisions to make about the ther-apeutic actions is left to the neurologists. Nevertheless, the signal-processing engineer should have some basics of neurology to be able to take a primary judgment about the goodness of his results.
The challenge in the field of neurology is the ignorance of the expected results. We cannot look inside the brain and compare our results to what we see. But we can assume that the results are as reliable as possible and help to make them even more reliable with today’s advanced signal processing algorithms and methods.
Chapter 3
Background
3.1
Magnetic Resonance Imaging (MRI)
3.1.1
Introduction
Magnetic resonance imaging (MRI) is an imaging technique based on the physical phe-nomenon of Nuclear Magnetic Resonance (NMR). It is used in medical settings to produce images of the inside of the human body. MRI can produce an image of the NMR signal in a thin slice through the human body. By scanning a set of such slices a volume of a part of the human body can be represented with MRI.
The principle of NMR is used in spectroscopy since the 1950s. The first human scan was taken in 1977 and the technology was adopted for clinical use around 1988.
3.1.2
Basic principles
NMR is the phenomenon that occurs when the sub-atomic particles in an atom nucleus, i.e. protons and neutrons, are immersed in a static magnetic field and exposed to a second oscillating magnetic field. If the net spin of all protons and neurons in an atomic nucleus is non-zero, the nucleus has a magnetic moment. This magnetic moment interacts with the magnetic fields in which it is immersed.
Hydrogen is the most abundant element in the human body (especially in water and fat). We consider a proton of a hydrogen atom to understand how particles with spin behave in a magnetic field. The magnetic moment of this proton causes the proton to behave like a tiny magnet with north and south poles. Placed in an external magnetic field B0
of a magnetic resonance scanner, the spin vector aligns itself with the external field like a magnet would. There are two possible energy configurations, a low and a high energy state where the spin vector is aligned parallel or anti-parallel with respect to the external magnetic field B0 as demonstrated by figure 3.2 (this figure as well as some others are
taken from Ola Friman’s dissertation [1]). Slightly more protons are in the low energy con-figuration, giving a net magnetic field in the same direction as the external magnetic field.
Figure 3.1: MR scanner (Philips)
Figure 3.2: The two different energy states a hydrogen spin vector can attain when placed in a static magnetic field B0
We have now spinning protons creating their own longitudinal magnetic field inside a high external magnetic field. In magnetic resonance, the interaction of the two magnetic fields results in precession of the protons.
The hydrogen protons are precessing with a frequency depending on the external magnetic field strength. This precessional frequency is
f = B0· γ (3.1)
B0 is the magnetic field strength measured in Tesla. γ = 42.58 MHz·T−1 is the
gyro-magnetic ratio of hydrogen. The precessional frequency f is also known as the Larmor frequency, after the man who discovered it.
If a radio frequency (RF) pulse of energy at the same frequency (the Larmor frequency of a hydrogen proton) is applied to these hydrogen protons for a fraction of a second, the energy will be transferred to the protons. The resultant magnetic vector points in a transverse direction. This process is called resonance and the protons are caused to precess in phase.
In the 1980’s the static magnetic field of most closed systems was ’only’ 0.5 T. The earth has a magnetic field strength of 30 to 60 µT. So the magnetic field of an 0.5 T MR scan-ner is approximately 10’000 times stronger than the earth’s magnetic field. In the 1990’s this was replaced by 1.5 T MR scanners and today 3 T systems are installed for medical use. Typical available MRI-Scanners have the following magnetic field strength inside the tube and the corresponding Larmor frequency of the second oscillating magnetic field (RF pulse magnetic field) to measure hydrogen atoms:
field [Tesla] 12 1 112 3 frequency [MHz] 21 43 64 128
In absence of this RF pulse the hydrogen protons gradually relax back, losing the trans-verse magnetism and their phase coherence. Individual protons relax at different rates, depending on what type of tissue the atoms are within. The precession of the protons causes an induction of a current in the receiver coils of the MR scanner. The RF pulses are repeated at a special interval (TR = time to repeat).
3.1.3
Imaging parameters
The electrical signals are analysed in terms of strength and location to create the MR image. The basic variations observed in the pulse sequences are T1 and T2 weighted
im-ages.
The time T1is the time of the longitudinal magnetization to recover, i.e. the rate at which
the longitudinal component of the magnetization vector regrows. Individual tissues have a unique T1-relaxation, some recover faster than others.
The time T2 is the time of the transverse magnetization to decay due to dephasing, i.e.
between the nuclei. T2-relaxation is also unique to individual tissues.
The nuclei in the studied ensemble are spatially distributed and may for this reason experience slightly different magnetic field strengths. The T2∗ time constant measures the combined effect of random nuclei interactions and magnetic field inhomogeneities. It naturally holds that T2 > T2∗. The current induced in the receiver coil due to the RF
pulse therefore mainly decays with a time constant T2∗. By applying different RF pulses in different sequences, all three time constants can be measured.
The resolution of the achieved images is normally 512 × 512 pixels or below. However as presented in the subsequent section, for fMRI the image size is normally 128 × 128 or even 64 × 64 pixels.
For more details the reader is directed to the large amount of literature on MRI (for example [9]).
3.1.4
Field of application
Radiologists (medical doctors specialized in the field of radiology) are trained in MRI to read the magnetic resonance images. MRI technologists operate the MRI scanner to obtain the images that the radiologist prescribes. Post processing algorithms are applied to magnetic resonance images to extract more information or enable better visualization of information in magnetic resonance images.
MRI has some advantages compared to other medical image registering techniques like Computer Tomography (CT) or other x-ray based imaging methods. First of all, the body is not exposed to x-rays that can harm the body. We just have to accept that persons with metallic implants that cannot be removed temporary like pace makers, tattoos with metal in the ink, permanent make up, metal plates, pins or other metallic implants are excluded of MRI examinations. For all other subjects, MRI does not cause any negative effects on the body. As far as we know it is a safe imaging method.
The disadvantages of MRI are for example the problem with claustrophobic people that should be scanned. Being in an MRI machine can be a very worrying situation for those people. The machine makes a tremendous amount of noise during the scan. Therefore patients are given earplugs or headphones during the examination. During the scans the patients are required to hold very still for a long period of time. MRI exams can range in length from 10 to 90 minutes. Involuntary motion artefacts like respiratory, cardiac and peristaltic movements can be a problem. Orthopedic hardware like screws can also cause artefacts on the images. MRI is a quite expensive technique.
MRI is generally used to visualize the soft tissues of the body. It has larger differences between soft tissues compared with CT. A typical field of application is MR angiography to observe blood vessels. Among many other applications is MRI ideal to diagnose mul-tiple sclerosis (MS), tumours, infections of the brain, spine or joints. It can also visualize torn ligaments in the wrist, knee and ankle or shoulder injuries.
3.2
Functional Magnetic Resonance Imaging (fMRI)
3.2.1
Introduction
Functional Magnetic Resonance Imaging (fMRI) is a technique to locate regions of the brain activated by a physical sensation or stimulus. This technique is harmless and pain-less and fMRI examinations can be carried out using advanced clinical MRI scanners available at all modern hospitals. For this reason this method of detection of brain ac-tivity attracts a growing interest in all fields of neurological research.
During an fMRI experiment, the brain of the subject is scanned repeatedly, using the fast imaging technique of echo planar imaging (EPI). Thus to perform an fMRI experiment an MR scanner with the ability to acquire EPI images is required. The T2∗ time constant is generally not utilized for the generation of conventional anatomical MRI images since it is susceptible to local magnetic field inhomogeneities. However, the T2∗ value is inter-esting in fMRI since it is related to neural activity.
The MR scanner repeats the RF pulse in an interval TR (time to repeat), during a period of 5 to 10 minutes. The in-plane size of the EPI images is 64 × 64 or 128 × 128 voxels and a stack of 10-40 slices are generally acquired. Around 100-200 such image volumes are repeatedly collected during the examination. The subject is required to carry out some task consisting of periods of activity and periods of rest while the images are acquired. The indicators of increased brain activity is the blood flow and blood oxygen concentration in the active area. Signal processing is used to reveal these regions of functional MRI scans. For better understanding, the physiology of the local blood flow changes in active brain areas are described below.
3.2.2
The Blood Oxygen Level Dependent (BOLD) response
We cannot observe directly changes in brain activation due to a stimulus. But already in the end of the 19th century were local blood flow changes in active brain areas pre-dicted. All neurons in the brain consume oxygen. The haemoglobin molecules in the blood provide the neurons continuously with new oxygen. An increase of neuronal ac-tivity increases also the demand of oxygen. This leads to an increased concentration of oxygenated blood in the capillaries surrounding the active brain area.
These local increases in blood flow can be mapped in fMRI as a change in raw image intensity. Oxygenated blood has different magnetic properties than deoxygenated blood. The T2∗ time constant becomes shorter in areas with low oxygen concentration and longer in areas with high oxygen concentration. In the active brain area the intensity of a voxel in the fMRI image is brighter (longer T2∗). This effect is called blood oxygen level dependent (BOLD) signal. This change of intensity can however not be detected by the bare eye since it is very small. Special signal processing techniques are needed.
Figure 3.3: Neural activity increases the blood flow in the active region to provide the neurons with more oxygen.
3.2.3
Brain activity detection
In a typical experiment of brain activity detection the subject in the MR scanner is ex-posed to a particular form of simulation. This simulation can be pictures shown by special glasses he wears, sound played by headphones or the movement of the subject’s finger. A series of low resolution MR scans is taken over time. The period between successive scans depends on the scanner’s performance. For some of the scans, the stimulus will be presented, for some other scans the stimulus will be absent. This activity/rest block is repeated throughout the whole experiment. There will be a BOLD response in the acti-vated brain areas due to the BOLD signal effect. To see which brain areas were actiacti-vated by the stimulus, the low resolution brain images are compared. The voxels whose time series contain a BOLD response have to be found to create a map of active brain areas.
3.3
Statistical Parametric Mapping (SPM)
3.3.1
Introduction
The map of similarity measures between a reference time series (the stimulus) and the voxel time series (fMRI data) is commonly referred to as a statistical parameter map. This construction of a spatially extended statistical map has been instantiated in a soft-ware called SPM (Statistical Parametric Mapping). SPM is an academic softsoft-ware toolkit for the analysis of functional imaging data [6].
The user should be familiar with the statistical, mathematical and image processing concepts of SPM. The SPM suite and associated theory was originally developed by Karl Friston to analyse PET data and made available in 1991 for the functional imaging community. By the time of this project, the current release is SPM2 (released the 12 of May 2003). It is designed to analyse data from fMRI, PET, SPECT and similar modalities.
3.3.2
SPM approach
SPM is a software consisting of Matlab functions, scripts and data files and some ex-ternally compiled C subroutines. One part of SPM is written to organize and interpret fMRI images. The image sequences can be realigned, spatially normalized into a
stan-dard space and smoothed by the software. In release SPM2, the statistical parametric mapping approach is voxel based, i.e. the parametric statistical model is assumed at each voxel.
After a spatial low-pass filtering, temporal correlation between a voxel time series and a hemodynamic model function is applied to each voxel to build a statistic image with maximum intensity projection (MIP). General Linear Model (GLM) parameters are used to test hypotheses (more about GLM see section 4.1).
3.3.3
SPM image format
SPM can deal with two or three dimensional data of any size (either image dimensions or voxel size). The data should be organized with a separate file for each scan.
SPM uses the header and flat binary image file format of ANALYZE [7]. A description of the main properties of the header and the image file is given below.
The image file
*.img: An uninterrupted array of (unsigned integer, signed short, signed integer, float or double) voxel values. Each *.img file has an associated header file that contains informa-tion about the image process.
To spatially normalize the images with SPM, we need to know the initial orientation of the images. The global variable sptl Ornt contains this orientation. After normalization with SPM the images have the following orientation:
X increases from Left to Right.
Y increases from Posterior to Anterior. Z increases from Inferior to Superior. The header file
*.hdr: The format of the 348 byte header file is that adopted by ANALYZE [8]. The necessary fields in the context of SPM include:
Field [SPM default global variable] image size {in voxels for x, y and z} [DIM] voxel size {in mm for x, y and z} [VOX] data type {see spm type.m [6]} [TYPE] a scaling coefficient {applied during memory mapping} [SCALE] offset of voxel values in *.img {in bytes} [OFFSET] the origin {(x,y,z) in voxels} [ORIGIN] description {a short string} [DESCRIP] It is important that these header files are correct. A common problem with using SPM usually reduce to incomplete or incorrect header files.
Chapter 4
Theory
To find neural brain activity we not only need background information about fMRI but also some signal processing theory. The fMRI process provides us four-dimensional data that we have to transform and adjust to see and interpret what we are searching for. Methods like the general linear model (GLM) detection are well known and applied in the problem of finding neural brain activity out of fMRI image series.
More detailed information is given if canonical correlation analysis (CCA) is used instead of GLM. If we want to evaluate the detection methods and their performance, statistical methods are needed. If we work with synthetic data to evaluate the methods, receiver operating characteristic (ROC) statistics may be applied for comparison. For real data sets those statistics cannot be applied. The location of the activity is not known in advance in real human data. The ROC has to be adapted. Therefore an introduction to modified ROC methods is given in this chapter.
4.1
GLM method
As mentioned above, the most commonly used neural activity detection method for fMRI is the General Linear Modelling (GLM). The idea is to set up a model for the stimulus applied to the subject in the MRI scanner and fit this model to the data. A good fit between this model (what we expect to see in the data) and the data itself means that the data was probably caused by the stimulation.
GLM is univariate, this means that the model is fit to each voxel time series separately. We only consider one voxel at the time. A univariate linear modelling is
y(t) = aT· x(t) + b + e(t). (4.1) y(t) is a scalar of the intensity value at each time point. x(t) is the model vector, also a function of time like y(t). a are the parameter estimates for x(t). b is a constant corresponding to the offset of the data. b can be easily included in a. e(t) is the error in the model fitting.
An estimation of the ”goodness of fit” for each voxel time series is then calculated. Sta-tistical tests on this estimation decide whether the voxel was affected by the stimulation or not. Or in other words; whether we find activation in the voxel or not.
4.2
Canonical Correlation Analysis (CCA)
In this section, an overview of Canonical Correlation Analysis (CCA) and constrained CCA are provided. The goal of this analysis is to find a relationship between two sets of variables. As the name implies, CCA uses correlation coefficients to quantify the relationship between the two sets of variables. The term ”canonical” pertains to the coordinate system in which the correlation is measured.
4.2.1
The canonical correlation
If we have two random variables x and y, the correlation coefficient ρ is a scalar between -1 and 1 measuring the degree of linear dependence between x and y.
ρ =def Cov[x, y]
pV [x]V [y]. (4.2) A strong dependence results in a correlation coefficient close to 1, ρ = 0 represents linear independence of x and y. For zero mean variables the relation
Cov[x, y] = E{[x − µx][y − µy]} (4.3)
reduces to Cov[x, y] = E[xy] (4.4) and V [x] = 1 n X (xi− µx)2 (4.5) to V [x] = 1 n X x2 = E[x2]. (4.6) Hence equation (4.2) is reduced to
ρ = E[xy]
pE[x2]E[y2]. (4.7)
In a CCA problem we have two multivariate variables, denoted x = [x1, . . . , xm]T and
y = [y1, . . . , yn]T. The problem consists of making linear combinations of the variables in
x and y. If either x or y is one-dimensional we essentially have the GLM approach. In a preprocessing step we have to remove the mean to have the condition of zero mean multivariate variables. Than we form two new scalar random variables x and y as linear combinations of the components in x and y,
x = wx1x1+ . . . + wxmxm = w T xx, y = wy1y1+ . . . + wynyn = w T yy. (4.8)
The goal of canonical correlation is now to find the linear combination weights wx =
[wx1, . . . , wxm]
T and w
y = [wy1, . . . , wyn]
T that give maximum correlation between x and
y. These weights are denoted regression weights. The equations (4.8) can be inserted in the definition of the correlation coefficient in equation (4.7),
ρ = E[(w
T
xx)(wTyy)]
q E[(wT
xx)(wxTx)]E[(wyTy)(wTyy)]
= w T yE[xyT]wy q (wT xE[xxT]wx)(wyTE[yyT]wy) = w T xCxywy q (wT xCxxwx)(wTyCyywy) , (4.9) Cxy is an (m × n) covariance matrix with the components of covariances between the
variables in x and y. Cxx and Cyy are the corresponding matrices for the components
of covariance between the variables within x and y respectively. CCA maximizes the expression (4.9) with respect to wx and wy. Hence the derivatives of equation (4.9) with
respect to wx and wy are set to zero, resulting in the following equation system,
Cxywˆy = ρλxCxxwˆx
Cyxwˆx = ρλyCyywˆy. (4.10)
λx and λy are scaling factors,
λx = λ−1y = s ˆ wT yCyywˆy ˆ wT xCxxwˆx . (4.11) They have no importance in the current context. We substitute ˆwx for ˆwy and vice versa
in equation (4.10) and obtain the following two eigenvalue problems, C−1xxCxyC−1yyCyxwˆx= ρ2wˆx,
Cyy−1CyxC−1xxCxywˆy = ρ2wˆy. (4.12)
The matrices C−1xxCxyC−1yyCyx and C−1yyCyxC−1xxCxy share the same eigenvalues. The
pri-mary solution, i.e. the first pair of canonical variates, is given by the eigenvectors wx
and wy belonging to the largest eigenvalue. It is recommended to solve the first equation
of the problem (4.12) to obtain wx and the canonical correlations. Then wy is found
through the second equation of problem (4.10). This approach ensures that wx and wy
4.2.2
Constrained CCA
In a general CCA problem, the regression weights wx and wy would be allowed to adopt
both negative and positive values. In the fMRI analysis context however the regression weights can be constrained or restricted to plausible values due to prior knowledge. Bet-ter detection performance can be obtained. It exists a procedure for constrained CCA where weights are restricted to be non-negative (Das and Sen [3]).
If we take Cxy to be the convolution matrix, we note that
ρ2unconstrained ≥ ρ2 constrained ≥ max i,j [Cxy]i,j 2 . (4.13) If the unconstrained solution has only positive components in wx and wy, this is the
solution also to the constrained problem. At the other end, the smallest squared cor-relation we are guaranteed to obtain is given by the largest corcor-relation between any of the input variables in x and y, i.e. the largest element in Cxy. It is shown that the
constrained solution equals an unconstrained solution to a modified CCA problem where one or several variables in x and y have been excluded. We will test all possible deletions to find the global optimum. This constrained version of CCA will significantly improve detection performance in fMRI (see chapters 6 and 7).
4.3
Principal Component Analysis (PCA)
An fMRI examination results in a large four dimensional spatio-temporal data set. The interesting information is hidden in the data. We search for a BOLD response present in some voxels to indicate them as active. In addition we will also search for a spatial shape of the active region. For both the exact BOLD response function as well as and the spatial extent of the active region we have little prior information. To build reasonable prototypes of the temporal and the spatial model we can use PCA.
Principal Component Analysis (PCA) is able to identify a subspace with interesting sig-nals within a multidimensional signal space. If we look at a BOLD response model, it can be the case that two parameters influence the shape of this hemodynamic response model in similar ways. We can vary the parameter values within reasonable limits and realize a corresponding model shape to each parameter set. Like this, we obtain samples of realistic BOLD responses.
PCA finds the eigen-timecourses ek(t) that minimize the mean square error
E X
t
(h(t) − ˆh(t))2 (4.14) for a partial expansion to P components:
ˆ h(t) = m(t) + P X k=1 wkek(t), (4.15)
where m(t) is the average simulated response,
ek(t) models the most significant deviations from the mean response,
P is the number of eigen-timecourses.
A candidate basis functions for the linear subspace model of the hemodynamic response are the P eigen-timecourses ek(t), k = 1 . . . P together with the average timecourse m(t).
4.4
Impulse response
To model the BOLD response we can use an impulse model based on the impulse re-sponse. The basic theory about the impulse response is presented in this section.
The impulse response is the response of a linear system to a unit sample impulse as system’s input. One of the suggested impulse responses is the gamma probability density function with the time parameter and two additional shape parameters τ and δ:
g(t, τ, δ) =t τ
δ
· e−δ(t−τ )τ (4.16)
Another used impulse function with five parameters excluding the time parameter is the Difference of Gamma probability density functions:
g(t, c, τ1, τ2, δ1, δ2) = t τ1 δ1 · e−δ1(t−τ1)τ1 − c · t τ2 δ2 · e−δ2(t−τ2)τ2 (4.17)
The BOLD response model can be generated by convolving a square wave function (box-car model) with one of these impulse responses (see chapter Methods 5.1.4).
4.5
Spatial filter functions
Spatial filters are used to take the spatial extent of the active region into account. This filters create weighted averages over the voxel values in the neighbourhood of the analysed voxel time series.
4.5.1
Gaussian filters
For the canonical correlation analysis we will need some spatial filter theory. We will use Gaussian filters. They are presented in this section.
The standard deviation σ of a Gaussian distribution of a variable X is defined as: σ = Σ(X − µ)
2
where µ is the mean of all realizations of X and N is the number of such realizations. This variance can also be expressed in terms of the full width of half maximum (FWMH) of the Gaussian distribution function:
σ = v u u t F W HM 2 2 2 · ln 2 . (4.19) A two dimensional circular Gaussian function in x and y with σ = σx = σy has the
following equation:
f (x, y, σ) = e−12 x2+y2
σ2 . (4.20)
In three dimensions the spherical Gaussian function in x, y and z with σ = σx = σy = σz
has the following equation:
f (x, y, z, σ) = e−12( r σ) 2 , (4.21) with r =px2+ y2+ z2.
4.5.2
Steerable filters
To model spatial extents of active brain regions with different orientation we will use steerable filters. The procedure to design a set of two or three dimensional steerable filter functions is described in this section.
1. A Gaussian smoothing filter f (z) with a certain size is defined. This Gaussian filter can also be replaced by any wished low pass filter for smoothing the fMRI data. This is the original filter.
2. The filter kernel is partitioned into four parts, one isotropic central part giso(z) and
three oriented parts gi(z), i = 1, 2, 3. These functions are used for weighting the
original filter. The four functions sum to one in every point. A linear combination of these weight functions can form any orientated filter shape.
The directions ˆni of the oriented weight functions for 2D filters are:
ˆ n1 = 1 0 , ˆn2 = 1/2 √ 3/2 and ˆn3 = −1/2 √ 3/2 The oriented weight functions are then obtained as
gi(z) = 4 3 1 − giso(z) zT − ˆni kzk 2 −1 4 , i = 1 . . . 3. (4.22) A natural choice of the center isotropic weight giso(z) is a Gaussian shaped kernel with
a width equal to half the original f (z) filter size. The weighting of the original Gaussian low pass filter with the four weighting functions gives the spatial filter basis functions,
Figure 4.1: Construction of a set of 2D spatial filters
fiso(z) = giso(z)f (z)
fi(z) = gi(z)f (z), i = 1 . . . 3
We can also construct three dimensional steerable filters with a minimum number of six oriented weight function with the following directions:
ˆ n1 = a 0 b , nˆ2 = −a 0 b , ˆn3 = b a 0 , ˆ n4 = b −a 0 , nˆ5 = 0 b a , ˆn6 = 0 b −a , a = p 2 10 + 2√5 , b = 1 + √ 5 p 10 + 2√5
The procedure for generating the filter basis functions is identical to the 2D case. The oriented weight functions for the three dimensional case are:
gi(z) = 1 − giso(z) zT − ˆni kzk 2 −1 6 , i = 1 . . . 6. (4.23)
4.6
Statistics
As in many other post-processing algorithms to fMRI data, our goal is to bring out signal from noisy background. Statistical methods will declare a brain voxel to be active (signal) or inactive (resting state). No matter which method we use to measure the performance of such an algorithm it is important to assess it in terms of accuracy (few voxels falsely declared to be active, also called specificity) and efficiency (few voxels falsely declared inactive, also called sensitivity). There exists one useful and popular tool for testing the efficiency of different diagnostic algorithms in fMRI. This tool is the receiver oper-ating characteristic (ROC) method. The diagnostic algorithms are applied to synthetic pseudo-human fMRI data and the ROC curve is produced to measure the efficiency of these algorithms.
In real human data the actual distribution of active voxels is not known. The ROC method has therefore to be modified to deal with real human data to be able to compare the efficiency of the fMRI detection methods like CCA.
In this section, ROC and the modified ROC method are presented.
4.6.1
ROC (synthetic datasets)
The task of fMRI is to designate a brain voxel as active or inactive based on the signal response at the voxel induced by a stimulus. Each voxel can take one of the following four conditions through a statistical test:
1. voxel truly declared to be active (Y ∩ T ) 2. voxel falsely declared to be active (Y ∩ F ) 3. voxel truly declared to be inactive (N ∩ F ) 4. voxel falsely declared to be inactive (N ∩ T )
Here T and F, respectively, stand for the known conditions of the voxel being signal (T) or resting data/noise (F), whereas Y and N, respectively, stand for the conditions of the voxel being declared active (Y) or inactive (N) through the statistical test.
The decision about a voxel’s condition is based upon the value of the test statistics variable X. This variable X is based on the post-processing algorithm and a function of the signal response X = f (ρ). The ROC curve used to measure the efficiency of the algorithm is a plot of P (Y |T ) against P (Y |F ) for different values of the variable X. P (Y |T ) is the conditional probability of declaring a voxel as active given that the response contains true signal for the observed value of X. P (Y |F ) is the conditional probability of declaring a voxel as active given that the response contains resting data or noise for the observed value of X. The ROC curve runs from (0, 0) to (1, 1). To characterize the power of a detection method to discriminate between active and inactive voxels we use the area under the ROC curve. The bigger this area, the better is the post-processing detection algorithm performance.
4.6.2
Modified ROC (real human datasets)
If the condition of a voxel being signal or resting data is not known, the test statistics variable X cannot be tested like described in the upper section. The conventional ROC curve cannot be obtained. Especially in the context of real human fMRI data we rarely know which voxels are truly signal (T).
P (Y |T ) and P (Y |F ) are denoted as true-positive fraction (TPF) and false-positive frac-tion (FPF) respectively. In the simulated case these two fracfrac-tions are known or can be estimated. In the case of real human data sets it is possible to estimate FPF from real resting-state fMRI data. We have to register pure resting data of the subject in whom the fMRI data is acquired. We can assume that this data is not related to activity, i.e. all voxels are inactive. Thus FPF can be estimated for different values of the test statistic X. The fraction of the resting data declared active is called fraction of resting positives (FRP).
TPF cannot be estimated from real fMRI data since we do not know which voxels are truly signal. Therefore we estimate the fraction of all voxels detected to be active for different values of the test statistic. This will include both true signal and noise. We can thus obtain the modified ROC curve for real data by using the fraction of detected active voxels instead of the TPF. This fraction is called fraction of active positives (FAP). The modified ROC curve is the plot of of FAP against FRP (have a look at the figures 6.1, 6.2 or 6.6 in the Results chapter 6). The conventional and the modified ROC curves are related to each other. Since we apply ROC to synthetic human data and modified ROC to real human data this relationship will not be explained in detail here. For more details about this relationship see article [2].
Chapter 5
Methods
The aim of this project is to evaluate and implement neural brain activity detection methods.
To evaluate the detection methods we can apply them on either synthetic test data or on real human data. The advantage of test data is that we can decide the region, geometry, size and strength of the activity signal. There exist well developed and easy statistical methods to compare the different detection methods if we work with synthetic data. In practice however neural scientists have to work with real human data. Real noise level and interference can never be exactly reproduced with a synthetic signal.
The detection methods and algorithms are programmed in Matlab for the evaluation task. The data has to be transformed into an appropriate format.
A great job of evaluation of test data with Matlab functions has already be done [1]. An overview of that work and results on synthetic data will be given in this and the following chapter. In addition the methods used to evaluate the detection methods on real data, in particular the CCA method, are presented within this chapter and the next chapter.
SPM is a widely used Matlab-based software available for free. It uses GLM to detect neural brain activity. The idea is to implement the Matlab functions developed under the evaluation task for CCA into SPM. SPM-users should be able to choose between these two different detection methods.The data used under the Matlab evaluation has to be adapted to SPM-format and the existing SPM code will be modified.
5.1
Evaluation of neural activity detection methods
In the following sections the reader will find information about the human data sets used to evaluate the performance of the detection methods on real fMRI data.
5.1.1
Brain data acquisition
The fMRI experiments were performed on an MR scanner at the University Hospital in Link¨oping, Sweden. The used MR scanner is an GE Signa LX Echospeed Plus 1.5 Tesla. Experimental setup
In the fMRI experiment performed for this thesis 120 image volumes are acquired with a sampling period TR of 4 seconds. The in-plane size of the EPI image slices is 64 × 64 pixels and a stack of 32 slices is acquired to create a volume. The slice thickness is 3 mm, the acquired matrix has the dimensions 64 × 64 × 32 voxels, the dimension of each voxel is 3 × 3 × 3 mm. During the experiment the test subject is instructed to perform a finger motion task. He is told to move his right digit during 40 seconds (or 10 scans) and rest during the following 40 seconds (10 scans).
In addition a series of pure resting scans is registered to have null data. It is supposed that there is no activity during those scans and all the signals within this data correspond to noise signals.
One file per brain slice scan is registered in a separate file. This is the raw fMRI data available for the analysis.
5.1.2
Brain data handling
Regardless of which kind of data, synthetic or real data, pre-processed or not, we have to transform it into the adequate format to handle it. In this thesis, three main formats turn up:
• The raw data of the MR scanner. We have one data file per scanned brain slice and per time point. The data is in our case floating point data with big-endian byte ordering.
• The *.mat data file used by the Matlab code. This data is an 4-dimensional matrix with the dimensions x × y × z × t.
• The *.img/*.hdr file data used by SPM. SPM uses the header and flat binary image file format of ANALYZE [7]. There is one data file per time point, i.e. one volume scan in each file. In our case the data is written in 32 bit floating point format. In the code available on demand, the reader finds all the necessary Matlab functions to transform the raw data of the MR scanner into the *.mat format of Matlab and the *.img/*.hdr of SPM. To facilitate the task of creation *.img/*.hdr files, a slightly modified version of the toolbox mri toolbox [] is used to agree with the specifications of SPM.
5.1.3
Implementation of the CCA method
This work is strongly related to the dissertation of Ola Friman [1]. Matlab code based upon his dissertation was available. I adapted the code to the desired requirements. An implementation of the CCA method into the existing software package of statistical parametric mapping, called SPM, is one of the aims of this thesis. The SPM software is a suite of Matlab functions and subroutines with some externally compiled C routines [6]. Therefore the implementation of all functions in Matlab was retained to facilitate the subsequent implementation of these functions in SPM.
Canonical correlation coefficient analysis (CCA) was not yet implemented in SPM. It is important to understand how CCA is implemented in Matlab to be able to implement it later in SPM. I will cover the important points of the work of Ola Friman [1] in this section. There will be some repetitions, but they are important to be mentioned to keep the integrity of this thesis.
Reasonable BOLD response models are needed as well as spatial filters for the CCA detection method. The following sections shows how the temporal and spatial modelling are used to produce the correlation map.
5.1.4
Generation of a hemodynamic response model
The response of neurons in the brain due to a stimulus is not exactly known. We have to make some simplifications to be able to simulate a feasible hemodynamic response. It is desirable to use a linear model of the hemodynamic response. We construct this response with a set of temporal basis functions. We want to have a linear model that captures the important variations in the hemodynamic response with as few temporal basis functions as possible to limit the cost of calculations and the risk of false positives.
Boxcar model
The most trivial BOLD response model is a boxcar model. The boxcar model has only one temporal basis function. The parameters we need to know are the sampling frequency (in Hz) the number of points per interval of resting state and activity period, respectively as well as the number of such intervals of resting and work.
Impulse model
Figure 5.2: Convolutional response model
A more sophisticated model is the impulse model. Again, we have just one temporal basis function. We can model the process of neural activation as a convolutional hemodynamic response model. We model a impulse response whose parameters are adapted to the hemodynamic response (see chapter Theory 4.4). A convolution of the impulse response with the boxcar model gives us the impulse model.
PCA generated model
Figure 5.3: PCA generated response function using a difference of Gamma functions. The mean response and the first eigen-timecourse are displayed.
To create a more realistic BOLD response model, we generate several samples of impulse responses changing the parameters in a sophisticated way. For each parameter of this convoluting model there is a range of values resulting in physiologically realistic shapes of responses. If we take the Difference of Gamma probability density functions
g(t, c, τ1, τ2, δ1, δ2) = t τ1 δ1 · e−δ1(t−τ1)τ1 − c · t τ2 δ2 · e−δ2(t−τ2)τ2 (5.1)
we vary the parameters c, τ1, τ2, δ1 and δ2 randomly within these ranges and produce a
large number of plausible response shapes gi(t). The mean values of the parameters are
chosen as c = 0.125, τ1 = 5.5, τ2 = 2τ1+ 5, δ1 = 6, δ2 = 14.5 [4].
A principal component analysis (PCA) applied to this set of BOLD responses gives us a set of orthogonal temporal basis functions. With those basis functions the hemodynamic response function can be reconstructed accurately with as few basis functions as possible. It was shown that it is normally sufficient to include the first eigen-timecourse in the subspace model. This will result in two hemodynamic basis functions, a mean response y1(t) = m(t) and the first eigen-timecourse y2(t) = e1(t) [1].
y(t) = wy1y1(t) + wy2y2(t) (5.2)
is the linear model for the hemodynamic response. Constraints to the temporal modelling
For a temporal model with more than one basis function, an arbitrary choice of weights would be able to produce an unrealistic BOLD response (see the examples in the grey zone in figure 5.4). It is possible to constrain the weights to give realistic shapes. We consider the PCA generated response model with two basis functions y(t) = wy1y1(t) + wy2y2(t). If
y1(t) and y2(t) are normalized to unit energy, the mean response y1(t) is obviously more
important than the first eigen-timecourse y2(t). So wy1 must be significantly larger than
wy2.
We will take the set of linear models generated to perform a PCA and make a least squares fit to plot them in a subspace spanned by y1(t) and y2(t).
As expected, the points lie around the mean response y1(t). y2(t) takes care of the small
deviations. We can now make a change of variable trick so that the weights of a realistic response are all non-negative. The change of variables
˜
y1(t) = y1(t) + αy2(t)
˜
y2(t) = y1(t) − αy2(t)
defines the cone in figure 5.4. If we realize a response y(t) = wy1y˜1(t) + wy2y˜2(t) only
realistic shapes lying in the first quadrant are obtained by constraining to non-negative weights. The parameter α ≥ 0 determines the angle of the cone, or equivalently how much weight we accord the eigen-timecourse y2(t) in relation to the main model shape
Figure 5.4: a) shows different BOLD model shapes. Realistic shapes lie around the main basis function y1(t). The change of variable trick is applied. Now the weights can be
constrained to non-negative values giving realistic shapes lying in the first quadrant in b)
5.1.5
Spatial modelling
To improve the detection of active voxels, we can exploit the spatial dependence among voxels. In general we have no prior knowledge about the shape of the active region. But it is more probable that a voxel is active if its neighbour shows activation too. The active brain areas where a BOLD response can be observed have in general an extent of several millimeters, thus a few voxels. If we just consider the temporal behaviour of a voxel it is possible that isolated voxels are declared to be active because of noise interference or artefacts without presence of real activity at this voxel.
Traditionally a fixed Gaussian filter was used for smoothing the images in a pre-processing step. This corresponds in the temporal model to use one single basis function as BOLD response model. If we have no information about the shape of the active region, it is recommendable to use more than one basis function to be able to shape every plausible region.
We will show in this section some shapes and models of spatial basis functions. They are the spatial counterpart to the temporal basis functions and used in the adaptive spatial filtering of the fMRI images.
Single voxel model
The simplest model for a spatial area is a voxel analysed in isolation of its neighbourhood. This is the case if we assume that the active region is smaller than the size of a voxel. This model is the one used if no spatial pre-filtering is applied to the fMRI image. Symmetric 3 × 3 model
Figure 5.5: Symmetric 3 × 3 model
We build the five symmetric basis functions of figure 5.5 for a symmetric 3 × 3 model. A linear combination of this basis functions can produce local activity patterns which are symmetric around the centre voxel. The shapes produced by this combinations become quite rough.
Scale adaptive model
Observations show that the spatial extent is the most obvious shape difference between active brain areas. As for the temporal model, we can produce several Gaussian filters by varying the shape parameter σ and apply a PCA to find a subset of basis functions able to reproduce a active region. Figure 5.6 shows the mean basis filter function and the first eigen-value basis filter found by PCA. This is an isotropic model.
It is also possible to apply a three dimensional Gaussian filter. As in the two dimensional case, we will apply a PCA to find a subset of spatial basis functions which account for the variations around the average shape.
Figure 5.6: Scale adaptive model with a subset of two basis functions
Figure 5.7: 2D orientation adaptive model Orientation adaptive model
In a general case the active regions have arbitrary orientations and can thus be captured better if we can adapt our model shape to different orientations and not just to the spatial extent as in the scale adaptive model. For this reason we use steerable filters (see figure 5.7).
The basis functions b2, b3 and b4 can be linearly combined so that the result is a rotated
basis function. Added to the isotropic part b1 smooth anisotropic shapes can be created.
b(ξ1, xi2) = wx1b1(ξ1, xi2) + wx2b2(ξ1, xi2) + wx3b3(ξ1, xi2) + wx4b4(ξ1, xi2). (5.3)
We have the ability to form anisotropic shapes in different orientations by changing the weights of the basis functions. To construct corresponding three dimensional shapes in any orientation, seven basis functions are required.
The price to pay compared to the scale adaptive model is the larger number of basis functions.
Constraints to the spatial modelling
Just as in the temporal modelling case an arbitrary choice of weights wxi in the case of
linear subspace models can result in unrealistic activity patterns. We will therefore im-pose constraints on the weights wxi to improve the detection performance. For example,
the weight of the centre voxel basis function in the 3 × 3 model, the original Gaussian shape in the scale adaptive model or the central isotropic basis function in the orientation
adaptive model has to be larger than for the remaining basis functions.
Again, we will use the change of variable trick presented in section 5.1.4 to be able to use non-negativity constraints for emphasizing the importance of the main basis functions.
5.1.6
Generation of a correlation map
We have one or more temporal basis function y and one or more spatial filter basis func-tion. Now we need to put the pieces together. We want to examine a voxel time series for neural activity.
In a first iteration the voxel time series neighbourhood is filtered with the spatial filter basis function(s). We obtain time series x, one for each basic function. If we filter the neighbourhood in 2 dimensions and with an orientation adaptive model we obtain 4 new time series. A 3 dimensional filtering with orientation adaptive filters results in 7 time series.
In a second iteration we try to find the best linear combination of these time series. We can simply apply a constrained CCA to the time series x and the temporal basis func-tions y to find the best linear combination weights ˆwx and ˆwy corresponding to the best
matching spatial filter and hemodynamic response.
To gain computational speed and for a better utilization of the constraints we can apply a 2-step method proposed in paper IV by Ola Friman [1]. First CCA is applied on the oriented filters and the temporal basis function y using the change of variable trick on y. Then CCA is applied a second time to the resulting time series and the central isotropic filter and the temporal basis function y. To give more importance to the centre neighbourhood currently under analysis a change of variable trick is not only applied to the y variables but also to the x variables as follows:
˜
x1(t) = xcenter(t),
˜
x2(t) = xcenter(t) + xoriented(t).
˜
x1(t) and ˜x2(t) are combined in the second step with positive regression coefficients to
give reasonable adaptive filter shapes. The possibilities offered by the constrained CCA are exploited in a better way compared to the direct one step procedure and computa-tional speed is gained.
At the end the correlation coefficient ρ between x and y is calculated for every voxel time series as described in chapter Theory 4.2.1. The coefficients ρ for every voxel gives the correlation map of the constrained CCA technique.
5.1.7
Noise filtering
The BOLD response is not the only signal present in the voxel time series. It is hard to determine which kind of noise is overlaid to the signal. Drifts can simply be removed by detrending the signal. We remove a continuous, piecewise linear trend of every voxel time series after having applied the spatial basis filters.
Effects of heart beat or respiration movements are not taken into account. It would be difficult to filter and remove them. They can cover the same frequency band as the BOLD response. We consider them just as noise. The signals of the experiments carried out for this project have a frequency of 1.3 bpm (bpm = beats per minute). We have the following values:
Cardiac signals : 60-80 bpm Respiratory signals : 15-20 bpm
BOLD response : 1-2 bpm (data set of this project)
I.e. For the experiments carried out for this project there is no risk that the BOLD re-sponse signal takes the same frequency as respiration movements or even as heart beating. But in the future or for other projects, if TR can be reduced and thus the frequency of the BOLD response signal can rise, the risk of wrong signal detection should be kept in mind. It can never be excluded that the test person is involuntary moving other parts of the body or thinking. This causes activation of other regions than the executed experiment is supposed to activate. We have to admit that such involuntary activation may cause uncorrelated hemodynamic responses to the desired and tested response function.
5.1.8
Statistics
The test statistics variable X are the correlation coefficients ρ acquired by one of the post-processing algorithms, either the GLM method or the canonical CCA method. ROC for synthetic datasets
The conventional ROC method requires knowledge about the distribution of true signal and resting state or noise. In the case of synthetic pseudo-human fMRI data, the data set consists of active regions as well as regions with resting state data. The locations of the regions presenting active signal are known.
The conventional ROC method can be applied to measure the efficiency of applied mul-tivariate post-processing algorithms to synthetic data sets. We can estimate TPF and FPF by calculating the number of voxels detected to be active correctly and incorrectly, respectively. We can than plot TPF against FPF to obtain the ROC curve.
Modified ROC for real human data sets
If we have real human data, the conventional ROC method cannot be applied. The loca-tion of the active voxels is not known in advance. Thus we will apply the modified ROC
method if we deal with real human data sets.
To estimate the fraction of resting positives (FRP) a pure resting state data set of the subject is needed. The subject has to refrain from any activity during the registration of the resting-state fMRI data. The detection algorithm is then applied with the same para-meters to the resting-state data as well as to the real data set in which activity is expected. We will apply constrained CCA to the resting-state data and to the real human data coming from an experiment to produce the correlation map, i.e. calculate the correlation coefficients ρ. The modified ROC curve is built with this data. The efficiency of a detection method is measured by the area under the modified ROC curve. The bigger this area, the better the detection method.
5.1.9
Matlab GUI
To compare and evaluate the different options to calculate correlation maps with CCA, I developed a small Matlab GUI. In addition the creation of a Matlab GUI was a good exercise for the later implementation in SPM and understanding of Matlab GUI’s (see figure 5.8).
Figure 5.8: Main Window of the Matlab GUI
The GUI is written to for the comparison of the detection methods. The user selects the fMRI images to work with (in *.mat format, see Appendix B). Then he chooses the spatial filter to apply as well as the hemodynamic response function to fulfil the constrained CCA. The correlation map is calculated with CCA and saved. In a last step the results (correlation map) can be displayed and compared to other correlation maps with modified ROC. The analysis can also be made with GLM and the corresponding map is saved as well.
5.2
SPM implementation
One goal of this project was to implement the CCA methods into the statistical pa-rameter mapping (SPM) described in chapter 3. Since we try to implement new code and functionalities into existing and well working code the following working steps are necessary.
1. Read and understand how the code is built.
2. Learn how the existing functionalities work by testing some examples. 3. Adapt the image format to SPM image format.
4. Decide how to include the CCA methods. 5. Write the new code and new functions.
6. Test the new functionalities. (Chapter results 6)
5.2.1
Structure of SPM software
SPM code is mainly written as Matlab code. GUI interfaces are constructed in order to give a more user friendly application. There exist two main functionalities. A first one to analyse PET and SPECT images, a second one for fMRI time-series images. We limit our studies to the second part for analyzing fMRI time-series.
fMRI time-series
The user interface to analyse fMRI time-series consists of tree main windows.
• A SPM menu window with buttons to choose the operations to handle the data. • A SPM interactive window with input/output handling between the user and the
functions.
• A SPM graphics window to show the results of the applied functions to the fMRI data.
Architecture
All the generated and needed analysis variables and parameters are saved within a special structure. SPM2 sets up a single structure (SPM.mat) at the beginning of each analysis. As the analysis proceeds, sub-fields are filled in. The key fields important to our appli-cation are (for more details see [6]):
SPM.xY : Data structure
SPM.nscan : Vector of scans per session SPM.xBF : Basis function structure
To implement constrained CCA we need to build a new SPM structure SPM.xCCA : CCA design structure.
5.2.2
Testing SPM2 code
A ordinary treatment of data with SPM consists of the following three steps: 1. Spatial pre-processing
The data is registered, warped and smoothed. 2. Estimation of model parameters
The GLM model parameters are estimated. The user has different options how these parameters should be estimated.
3. SPM
Statistical parameter mapping. The maximum intensity projection (MIP) can be displayed and the models used can be visualized.
SPM is a quite advanced software, therefore a lot of options are available within all of this three steps. It is left to the user which of the options he needs for his analysis purpose and which are redundant.
5.2.3
Integration of CCA methods
To include the CCA methods into the existing SPM software area we have to determine some requirements.
• The CCA method is used in the context of fMRI analysis. It will be a part of the fMRI structure of SPM.
• The existing functionalities have to work as supplied before. No changes will be made to the existing functionalities.
• If one wishes to work with the CCA methods it should be easy to access to the new code. We will write a toolbox to overload some SPM functions and add new SPM functions performing the constrained CCA task.
• The architecture has to be conserved as well as possible. A lot of variables of SPM are not needed by the CCA method. The structure (SPM.mat) should still contain as many as possible of the variables in order to fulfil the ordinary SPM functions. SPM code compatible with those requirements was implemented. The functionalities of the new code are presented in chapter 6 and a guideline is given in the appendix.
Chapter 6
Results
In this chapter we will show how the methods described in chapter Methods 5 work.
6.1
Some important values
The fMRI data from the right hand finger tapping experiment with isotropic voxel size is used for all of the comparisons and tests in this chapter. Some important values and characteristics of the used fMRI are mentioned below:
TR = 4 [s]
data matrix = 64 × 64 × 32 [voxels] voxel size = 3 × 3 × 3 [mm]
t = 120 [scans] duration of onsets = 10 [scans]
vector of onsets = [6 26 46 66 86 106]
For the constrained CCA analysis in SPM a correlation threshold of 0.55 was used to display the correlation maps.
6.2
Evaluation on synthetic data
As mentioned in the Methods chapter 5, Ola Friman already evaluated the CCA method on synthetic data in his work and papers. In this work his basic results are used and completed with the results on real human data.
6.2.1
Related work
The GLM approach to detect neural brain activity is equivalent to an unconstrained CCA with just one spatial basis filter. By applying CCA, this approach is generalized and aug-mented with the use several spatial basis functions to improve the detection performance. The spatial filter basis functions are the spatial counterpart of the use of temporal basis functions for modeling the hemodynamic response.
It has been shown that the following models can improve the neural brain activity detec-tion results [1].
PCA generated BOLD response model
The better the model for a hemodynamic response function matches the the BOLD response function, so much the better is the efficiency to find activity in a voxel time-series. The commonly used model is a convolution between a binary reference function (boxcar model) and an impulse response consisting of a difference of two Gamma functions. If PCA is applied to a sample set of randomly generated responses with parameter values chosen within reasonable ranges it is possible to obtain well fitting linear model shapes. It has been shown that the first principal component accounts for about 80% of the variations, therefore a model consisting of two basis function, y1(t) the average response
and y2(t) the first principal eigen-response is adequate. If we make the change of variable
trick, the weights giving realistic response shapes are non-negative and can be found by constrained CCA.
Steerable spatial filters
Steerable filters can adapt to both orientation and size of an active brain area. Compar-isons on synthetic test data show that the detection of all kind of active region geometries are superior if we use steerable spatial filters instead of scale adaptive or other simple spatial filters.
Constrained analysis
The results of CCA using constrained and unconstrained regression weights show that the constrained analysis with adaptive filters gives more distinct correlation maps with better contrasts. Without constraints the number of active voxels detected is bigger due to the large freedom in constructing suitable filters. It was shown that the non-negativity constraints can suppress the correlation levels in non-active areas and improve specificity. 3D spatial filtering
3D filtering improves the possibilities to detect active areas. However the price to pay is an increase of computation time.
6.3
Evaluation on real human data
In the following sections the results of the same techniques applied to real human data are shown. The performance of the detection methods and techniques are measured using the modified ROC method. The correlation maps have been generated using SPM with the new CCA code added.
Figure 6.1: Modified ROC curve of scale versus orientation adaptive filters
6.3.1
Steerable spatial filters
We will show the advantage of using orientation adaptive filters compared to scale adap-tive or even fixed spatial filters. Steerable basis filters are used to adapt the size and orientation of the active brain area. By generating the modified ROC curve of figure 6.1 we see how the orientation adaptive filters are more specific in detecting active brain areas compared to an scale adaptive model.
6.3.2
Constrained CCA
Non-negativity constraints are added to the hemodynamic model as well as to the spatial filter model to avoid unrealistic model shapes. The α-value for the change of variable has been set to ±0.4. The constraints on the spatial filter model are automatically calculated by SPM.
As for the synthetic data the modified ROC curve of figure 6.2 clearly shows the improved specificity of detecting active voxels in fMRI data.
