A Matlab Toolbox for fMRI Data Analysis: Detection, Estimation and Brain Connectivity

(1)

Master thesis

A Matlab Toolbox for fMRI Data

Analysis: Detection, Estimation and

Brain Connectivity

by

Kiran Kumar Budde

LiTH-ISY-EX--12/4600--SE

05-09-2012

(2)

(3)

Master thesis

A Matlab Toolbox for fMRI Data

Analysis: Detection, Estimation and

Brain Connectivity

by

Kiran Kumar Budde

LiTH-ISY-EX--12/4600--SE

05-09-2012

(4)

(5)

Abstract

Functional Magnetic Resonance Imaging (fMRI) is one of the best techniques for neuroimaging and have revolutionized the way to understand the brain functions. It measures the changes in the blood oxygen level-dependent (BOLD) signal which is related to the neuronal activity. Complexity of the data, presence of different types of noises and the massive amount of data makes the fMRI data analysis a challenging one. It demands efficient signal processing and statistical analysis methods. The inference of the analysis are used by the physicians, neurologists and researchers for better under-standing of the brain functions.

The purpose of this study is to design a toolbox for fMRI data analysis. It includes methods to detect the brain activity maps, estimation of the hemo-dynamic response (HDR) and the connectivity of the brain structures. This toolbox provides methods for detection of activated brain regions measured with Bayesian estimator. Results are compared with the conventional meth-ods such as t-test, ordinary least squares (OLS) and weighted least squares (WLS). Brain activation and HDR are estimated with linear adaptive model and nonlinear method based on radial basis function (RBF) neural network. Nonlinear autoregressive with exogenous inputs (NARX) neural network is developed to model the dynamics of the fMRI data. This toolbox also provides methods to brain connectivity such as functional connectivity and effective connectivity. These methods are examined on simulated and real fMRI datasets.

(6)

(7)

Acknowledgement

I would like to express my deepest gratitude towards my supervisor Dr. Sadasivan Puthusserypady for giving me an opportunity to work in the field of functional magnetic resonance imaging. I would like to mention that his support is invaluable and have patiently led me throughout my work. With-out his constant support this thesis work would not be finished.

I would like to thank Dr. Luo Huaien for providing his algorithms for fMRI Matlab toolbox design and sharing his knowledge during the project. I thank my examiner Dr. Maria Magnsson and Dr. G¨oran Salerud for the thesis proposal and allowing me to do the thesis at the Technical Uni-versity of Denmark (DTU). I appreciate the help from Dr. G¨oran Salerud the starting of my master’s programme.

I would like to thank all those who provided invaluable advices and the help during my masters at Linkoping University and Technical University of Denmark.

Last but not least, I would like to say that the support from my family is immeasurable under all the circumstances.

(8)

Abbreviations

AIC Akaike Information Criterion

AR Autoregressive

ARX Autoregressive Model with Exogenous Inputs BOLD Blood Oxygenation Level Dependent

DCT Discrete Cosine Transform DWT Discrete Wavelet Transform

EEG Electroencephalography

EPI Echo Planar Imaging

fMRI Functional Magnetic Resonance Imaging

FPR False Positive Ratio

GLM General Linear Model

GLM Graphical User Interface

HDR Hemodynamic Response

HRF Hemodynamic Response Function

ISI Inter-Stimulus Intervals

LMS Least Mean Square

LS Least Squares

MLP Multi-Layer Perceptrons

MVAR Multivariate autoregressive

MR Magnetic Resonance

MRI Magnetic Resonance Imaging

NARX Nonlinear Autoregressive with Exogenous Inputs

OLS Ordinary Least Squares

PET Positron Emission Tomography

RBF Radial Basis Function

ROC Receiver Operator Characteristic

SNR Signal-to-Noise Ratio

SPM Statistical Parametric Mapping

TPR True Positive Ratio

(9)

Introduction

1.1 Introduction

Brain is the most fascinating, mysterious and least understood organ of the human body. For the last few years, functional brain imaging techniques have been advanced tremendously. For understanding the brain functions and brain mappings, a powerful tool like functional magnetic resonance (fMRI) can be used. fMRI measures the changes in blood oxygen level-dependent (BOLD) signals which are related to neural activity [1].

The purpose of this thesis is to develop a Matlab toolbox for the fMRI data analysis. It includes methods to detect the brain activation, estimation of hemodynamic response (HDR) and the connectivity of brain structures.

The major features of the toolbox are:

• To construct a flexible design matrix in the general linear model (GLM) under the Bayesian framework. The Bayesian approach is extended to nonstationary noise and drift models. These frameworks provide accurate detection and avoid multiple comparison problems in con-ventional methods. This estimator detect more real activation of sim-ulated and real fMRI datasets when compared with the traditional methods such as t-test, ordinary least squares (OLS) and weighted least squares (WLS).

• To provide methods for estimation of brain activation and HDR for linear and nonlinear properties of event related designs. The linear and non-linear properties are based on the inter-stimulus interval (ISI).

(12)

1.2. OUTLINE OF THE THESIS

• Nonlinear method based on radial basis function (RBF) neural network to detect the spatial activation.

• Nonlinear autoregressive with exogenous inputs (NARX) neural net-works to model the dynamics of fMRI data.

• It describes methods on functional connectivity analysis of brain re-gions using non-linear cross correlation analysis. Moreover, it mea-sures the directional interactions between spatially separated neural populations by using the Granger causality.

1.2 Outline of the thesis

Chapter 2: Briefly explains fMRI basics and discuss various aspects which are related to this project.

Chapter 3: This chapter describes the design of the toolbox including vi-sualization and loading of experimental real and simulated data. Different brain activity detection methods are explained in the presence of noise and drift.

Chapter 4: This chapter provides different nonlinear methods for detection of brain activation and estimation of the HDR.

Chapter 5: This chapter describes methods to investigate brain connec-tivity structures: functional connecconnec-tivity and effective connecconnec-tivity. The functional connectivity is investigated using the nonlinear cross correlation and effective connectivity is investigated using the Granger causality. Chapter 6: This chapter gives the conclusion of the thesis and the scope for future work.

(13)

Chapter 2

Background

2.1 Functional magnetic resonance imaging

The fMRI is one of the most advanced neuroimaging techniques which uses the standard magnetic resonance imaging (MRI) to examine the brain func-tions. It is a widely used technique because of its better spatial resolution when compared to electroencephalogram (EEG) and better temporal res-olution compared to positron emission tomography (PET). Moreover, it is non-invasive, gives non-ionizing radiation and has high sensitivity [2].

2.1.1 BOLD signal generation

As shown in Figure 2.1, when neural activity in the brain increases, neurons consume more oxygen and demand more oxygen. This results in increased blood flow at the activated areas. As a result, the oxygen concentration increases and decreases the deoxyhemoglobin. The changes in oxygen con-centration level alters the main magnetic field because the oxyhemoglobin is diamagnetic and deoxyhemoglobin is paramagnetic. T2*1_{time becomes}

shorter at low oxygen concentration and higher at high oxygen concentra-tion areas. Hence, the MR signal depends on oxygen level. This effect is referred as the Blood Oxygen Level Dependent (BOLD) effect [3].

1_{MR images contrast determined by properties of tissue being imaged, different tissues}

have different relaxation times T1,T2 and T2*. T2* relaxation time is time constant that describes the exponential decay of signal, due to spin-spin interactions, magnetic field inhomogeneities and susceptibility effects. T2* weighted imaging is commonly used in fMRI [4].

(14)

2.1. FUNCTIONAL MAGNETIC RESONANCE IMAGING

Figure 2.1: Block diagram of BOLD signal generation.

2.1.2 Hemodynamic response (HDR)

The change in MR signal on T2* images triggered by the neuronal activity is known as the HDR. It can result from the reduction in the amount of de-oxygenated blood and it represents temporal properties of brain. Figure 2.2 shows the sketch of a typical HDR. Its shape varies with activation: ampli-tude increases with rate of neural activity and width increases with increase in duration of the neuronal activity [5]. The HDR can be summarized in three phases [6].

(15)

Figure 2.2: Schematic representation of HDR.

• Intial dip: During the initial short time (1-2 sec), the MR signal de-creases below baseline after beginning of neural activity. It is caused by the transient increase in deoxyhemoglobin due to the oxygen con-sumption.

• Overcompensation: After the initial dip, more oxygenated blood is supplied to the area than extracted and deoxygenated hemoglobin decreases due to the increase of neuronal activity and results in MR signal increase above baseline at about 2 to 5 seconds.

• Undershoot: After finishing neuronal activity, the MR signal am-plitude gradually decreases below the baseline level and reaches the baseline level due to a combination of decreased blood flow and in-creased blood volume.

2.1.3 fMRI Experimental design

There are two schemes of experimental designs that are generally used for fMRI experiments: block design and event-related design [7].

(16)

design for detection of brain activity. It provides large signal to noise ra-tio (SNR). Figure 2.3 represents a schematic diagram of the BOLD signal (blue color) and the block design experiment (magenta color). The square waveform represent active (stimulus ON) and rest (stimulus OFF) of brain activation. This is not the best design for temporal activity estimation. In event related design, discrete and short duration events are presented one at a time and separated by random ISI that can range depending on the exper-iment. Figure 2.4 shows event related BOLD signal with the corresponding timing of events.

(17)

Figure 2.4: Event-related design BOLD signal.

2.1.4 fMRI experimental data used in this thesis

In this thesis, two types of experiments are examined. One is the block design (DATA BLOCK) and another is the event related design (DATA EVENT). The block design dataset was obtained from the Statistical Para-metric Mapping (SPM) website [8].

DATA BLOCK

The real fMRI experiment was designed for the activation of the auditory areas. The functional data or Echo Planar imaging (EPI)images acquired using 2T Siemens MAGNETOM Vision system. The experiment data set and its details are available on [8]. During the experiment, bisyllabic words are presented binaurally. The first few scans are discarded due to T1 effects then starts from fm00223 004 image. The total number of acquisitions are 96 and it is divided into 16 blocks and each block contains 6 acquisitions. In this experiment alternate condition of the baseline (rest) and activation are applied for an extended period of time.

(18)

2.2. FMRI DATA ANALYSIS

2.2 fMRI data analysis

In fMRI scanning, the whole brain or part of the brain is scanned over time period and generates 4D data or sequence of 3D MR images.

2.2.1 Preprocessing

The purpose of preprocessing is to remove the unwanted data to prepare for the statistical analysis and enhance the brain activity mappings. There are several steps to be followed [9]. In this thesis the following preprocessing steps are performed by the SPM software [10]:

• Realignment: During the data acquisition, intensity of the voxel of resting and activation assumes that no motion of subject has occurred. The motion is caused by the subject movement due to long period of scanning. This problem can be reduced by the realignment procedure with estimation and correction [11].

• Co-registration: The low resolution fMRI images are aligned with the anatomical MRI images, which can be images of the same subject or a standard template [12].

• Normalize: It is a procedure to wrap the functional data onto a coordinate system or template space [13].

• Smoothing: The fMRI images are smoothed across adjacent voxels to improve the results of brain mapping and to increase the SNR [14].

2.3 fMRI data modeling

2.3.1 Temporal modeling

From an engineering point of view, the fMRI data analysis can be considered as the analysis of the response of the system to a given input. As shown in Figure 2.5, the system is the human brain and measuring device, input is the experimental design and output is the observed BOLD signal.

(19)

2.3. FMRI DATA MODELING

To understand the complexity of fMRI data, spatial and temporal prop-erties are required. The measured fMRI signal or voxel time series is the combination of BOLD signal, drift and noise, i.e. the change in voxel inten-sity represent whether the BOLD signal is active or not.

Measured fMRI time series = BOLD signal + drift + noise. (2.1)

2.3.2 BOLD Model

This model assumes a linear system, i.e. the BOLD response expresses the convolution of the stimulus function and impulse response [15],

yb(t) = h(t) ⊗ s(t) =

Z ∞

0

h(τ )s(t − τ )dτ, (2.2) where ⊗ denotes the convolution operation, s(t) is the stimulus function and h(t) is impulse response function and it is called hemodynamic response function (HRF). The HRF waveform is modelled by poisson, Gaussian func-tion and difference of gamma funcfunc-tions [16]. Here, the widely used difference of gamma function model is used and it can be represented by,

h(t) = t d1 a1 exp −(t − d1) b1 − c t d2 a2 exp −(t − d2) b2 , (2.3) where di = aibirepresents the amplitude of the peak and c represents the

un-dershoot. The normally used parameters are a1= 6, a2= 12, b1= b2= 0.9,

and c = 0.35 [17].

2.3.3 Noise and Drift

The BOLD signal is influenced by strong random noise which results in low SNR. The sources of noise are scanner induced noise and subject movements during scanning such as head and lower jaw movements. The residual of the imperfect model is also one of the noise components [15].

Drift is another component that disturbs the fMRI time series. It shows up as a slow varying interference or trend in the fMRI time series. This drift often comes from physiological processes like the cardiac and inspira-tion process [18] as well as from instability of the magnetic field [19]. The common way to remove the drift in fMRI time series is either by using a high-pass filter or a drift model.

(20)

Chapter 3

Detection

This fMRI Matlab toolbox provides different methods for detection of ac-tivated brain regions and estimation of the HDR. The proposed toolbox consists of four functionalities. The first section describes how to display the cross sectional brain images and to load the experimental data. The sec-ond section provides different brain detection methods. The third section provides methods for estimation of the HDR and the fourth section describes brain connectivity algorithms. Figure 3.1 shows the graphical user interface (GUI) for the designed toolbox.

The first section explains the display of MRI cross sectional images and loading of real and simulated fMRI datasets. Second section explains de-termination of flexible design in the GLM through the Bayesian learning procedure and it provides accurate activation compared with conventional t-test method. The Bayesian estimator is extended to non-stationary noise model and it provides accurate activation results when compared to OLS and WLS. Moreover, the Bayesian estimator can be extended to remove drift component present in the fMRI time series. Studies on simulated and real fMRI data show that the Bayesian estimator accurately detects the brain activated regions to specific task. The detection and estimation algo-rithms developed by Dr. Luo Huaien have been adopted in the design of this toolbox.

(21)

Figure 3.1: This toolbox GUI comprises 3 windows, i) Main window, ii) Subwindow, iii) Graphical window. The main window consists four sections, i) Processing and visualization, ii) Detection of activated brain regions, iii) Estimation of HDR and iv) Analysis of brain connectivity.

(22)

3.1. LOADING AND VISUALIZATION

3.1 Loading and visualization

3.1.1 Display

The fMRI datasets are provided in the analyze file format. It is an image file format for storing MRI data and consists of two types of files, an image file and a header file, with .img and .hdr extensions, respectively. An image file contains uncompressed pixel data or a set of cross sectional images. The header file contains history and dimensions of the data. Before the fMRI data analysis starts, the dataset is preprocessed by the SPM software to apply realignment, coregistraton, normalization and smoothing. In Figure 3.2, the graphical window displays the brain cross-sectional images and the subwindow provides information about the images.

(23)

3.1.2 Loading real fMRI dataset

To load real fMRI data set, the toolbox requires series of MRI scanned data or session data1_{, which are stored in the analyze data format (hdr and img).}

Figure 3.3 shows the result of loaded auditory fMRI activation data. It also provides dataset information such as the number of slices, image and voxel dimensions. The graphical window shows only one image and other data details. The loaded image files in Figure 3.3 are from

swrfM00223 016.img, swrfM00223 016.hdr upto swrfM00223 099.img, swrfM00223 099.hdr.

Figure 3.3: Loading real fMRI dataset.

3.1.3 Loading simulated fMRI dataset

To load simulated data, simulate the activated regions at particular posi-tions of one image slice, by adding the BOLD signal repeatedly over time. As shown in Figure 3.4, the BOLD signal (block design or event related) is added at particular positions in the image. It is also possible to simu-late fMRI data with noise, such as time varying and fractional noise. In this thesis, two types of datasets are simulated: DATA BLOCK and DATA EVENT. DATA BLOCK is simulated with block design signal and DATA EVENT is simulated with event related signal.

(24)

In GUI, select the Load Simulated button. It opens a subwindow. Then choose the type of design (DATA BLOCK or DATA EVENT), type of noise (time varying or fractional) and then simulate it. In the graphical window, simulated data with activated positions are displayed (Figure 3.4) and noisy images are displayed in Figure 3.5.

(25)

Figure 3.5: Simulated image with time varying noise.

3.1.4 Loading stimulus

To load the stimuli or experimental paradigm, we need to have knowledge about the experiment. For block or event-related design experiments, several values must be entered: active and rest block length, number of blocks, repetition time (TR), signal duration and duration of rest. Block design experiment data and information are available in [4]. Figure 3.6 shows the block design paradigm of auditory experiment and Figure 3.7 shows the event-related paradigm.

(26)

Figure 3.6: Block design paradigm.

(27)

3.2. DETECTION OF ACTIVATED BRAIN REGIONS

3.2 Detection of activated brain regions

3.2.1 Flexible design matrix

This section briefly explains the general linear model (GLM). Moreover, the construction of a design matrix in GLM using Bayesian analysis is considered and compared with the conventional methods like the t-test [20].

3.2.1.1 t-test

The t-test is a conventional test used in brain mapping. It compares the av-erage active condition with the avav-erage rest condition. The t-test is defined as [21] t = _qX1− X2 S2 1 n1 + S2 2 n2 , (3.1)

where X1, S12 and n1 are the mean, variance and sample number for the

active fMRI samples. For the resting fMRI time samples X2, S22and n2are

the mean, variance and sample number.

Figure 3.8: Detection results of simulated fMRI data using the t-test method.

(28)

Figure 3.9: Detection results of real fMRI data using the t-test method. Figure 3.8 and 3.9 show the detection results of the t-test method applied to the simulated and real fMRI datasets. It can be clearly observed from that the t-test method also detects false detection. Therefore for the real fMRI data, analyzing the performance with this method makes it difficult. This is because we lack a reference which could serve as the true activation of the brain.

(29)

3.2.1.2 Flexible design matrix with sparse Bayesian method General linear model (GLM)

The GLM is a fundamental method for fMRI data analysis. In this, the measured fMRI time series is described as the combination of regressors and their vector parameters [22].

y = Φw + , (3.2)

where y is an N × 1 measured fMRI time series, Φ is an N × P ma-trix known as the design mama-trix, w is a P × 1 vector of estimated unknown parameters and is N × 1 noise vector. The design matrix contains all avail-able knowledge about the experiment. In the design matrix each row (N ) represents one time point (one scan) of the regressor and column (P ) repre-sents the corresponding regressor or explanatory variable like the canonical BOLD response, constant one vector and discrete cosine transform (DCT) basis functions [23].

Sparse Bayesian Learning Algorithm

In the classical approach, the design matrix is not flexible which may lead to problems. More regressors leads to model over-fitting problems and few basis functions cannot filter out the noise. The selection of regressors in the design matrix is a critical issue. The Bayesian approach is used to construct a flexible design matrix. During the learning procedure, the data itself es-timate the regressors. It determines which regressors are useful and ignores the ones that are irrelevant [20].

This algorithm is developed based on the assumption that the noise is stationary. Initially, the design matrix is initialized with the BOLD signal, a constant vector of value 1 and radial basis functions. Throughout the learn-ing procedure, the regressors and their weight coefficients are estimated and irrelevant regressors are discarded.

Figures 3.10 and 3.11 show that the Bayesian algorithm provides better results for both simulated and real fMRI datasets. The performance is better than the t-test method. Both the t-test and Bayesian estimator results clearly show that auditory activated areas are identified. As mentioned before, the sparse Bayesian algorithm assumes that the noise is stationary. However, the fMRI data also contains nonstationary noise which could be removed by nonstationary noise models.

(30)

Figure 3.10: Detection results of simulated fMRI data using the sparse Bayesian learning method.

Figure 3.11: Detection results of real fMRI data using the sparse Bayesian learning method.

(31)

By observing the results on real fMRI data, it is clearly seen that the au-ditory cortical areas has been identified. The Bayesian estimator detects the more real activated brain regions compared to conventional t-test method.

Figure 3.12: ROC curve of simulated data of the t-test and the Bayesian estimator.

Receiver operator characteristic (ROC) analysis is used to investigate clear comparison of detection ability of the t-test and the Bayesian estimator methods. The ROC curve plots the true positive ratio (TPR) against the false positive ratio (FPR) [24]. From Figure 3.12, we can say that, the Bayesian estimator is better compared to t-test.

(32)

3.2.2 Nonstationary noise models

The noise sources are scanner induced noise, physical movements of the sub-ject such as lower jaw movements, and neurophysiological processes such as number of neurons involved in specific tasks at different time points and the background memory process. The noise sources may cause the noise variance to be changed [25]. The aim of fMRI data analysis is to distinguish the activated BOLD signal from noisy data and to determine the activated regions for particular tasks. The high amount of noise in fMRI makes it difficult and challenging for data analysis.

This section briefly explains the Bayesian estimator for detection of acti-vated regions in the case of time varying noise. The performance is compared with the OLS estimator and the WLS estimator [26].

3.2.2.1 OLS estimator

As introduced in Section 3.1, the General linear model is defined as

y = Φw + , (3.3)

where ynand nare the observed fMRI time series and noise at the nthvoxel

respectively. n is assumed to be independent and identically distributed

(i.i.d) white noise. The least square estimator of the parameter w is defined as

ˆ

wn(OLS)= (ΦTΦ)−1ΦTyn. (3.4)

The assumption of i.i.d white noise is an inappropriate assumption re-garding the nonstationary nature of fMRI signals. The design matrix should be a full rank square matrix. If the design matrix is a rectangular matrix, the matrix inversion is not possible. Therefore in order to find the matrix inversion, firstly we multiply the design matrix with its transpose, which results in a square matrix [22]. The activated voxels are calculated by using static student distribution t-test

t = c T_w_ˆ p cT_Λ ˆ wc , (3.5)

where c is the contrast vector and Λwˆ is the covariance matrix of the

esti-mated parameter ˆw.

The Figure 3.13 shows the detection results of a simulated fMRI data using the OLS method in presence of time varying noise, here, the noise range is from -0.4 to 6.7 dB. The OLS method detects the activated regions as well as false detection. From Figure 3.14, it is clearly shows the detection of auditory cortical areas by using this method.

(33)

Figure 3.13: Detection results of simulated fMRI data using the OLS esti-mator.

(34)

3.2.2.2 WLS estimator

The noise variance of time varying noise model is stated to increases in a multiplicative fashion [27]. For this reason the noise variance is modeled with a scaling matrix and overall variance. The noise model is assumed to be time dependent, with precision matrix Bn. The precision matrix is the

inverse of the covariance matrix and it is defined as [25],

Bn =      s1 0 . . . 0 0 s2 . . . 0 .. . ... . .. ... 0 0 . . . sT      βn= Sβn, (3.6)

where s1, s2, ..sT are the scaling parameters and T is the number of samples

in the fMRI time series. Here, S is the scaling diagonal matrix and Bn is

the noise precision at nth_voxel.

The WLS estimator of parameter w is defined as ˆ

wn(WLS)= (ΦTSΦ)−1ΦTSyn. (3.7)

The WLS estimator requires accurate estimation of the scaling matrix. In traditional methods, the residual of the OLS estimator is used. The residual of the OLS estimator is given by

rn= yn− Φ ˆwn(OLS)= yn− Φ(ΦTΦ)−1ΦTyn. (3.8)

The overall precision of the nth _{voxel time series and the inverse scaling}

matrix are calculated as in [17], ˆ βn= T − rank(Φ) rT nrn , (3.9) ˆsinv= PN n=1diag(βnrnrTn) N , (3.10)

where, N is the total number of voxels and the operator diag(.) creates the diagonal of a square matrix into a column vector.

This method is simple to implement. The residual is equal to Sβn and

then is considered as unbiased estimator. The WLS method uses the residu-als of the OLS method as a covariance matrix to estimate the scaling matrix. It fails to estimate an accurate scaling matrix. The time varying Bayesian estimator, however, is able to accurately find the scaling matrix.

(35)

Figure 3.15: Detection results of simulated fMRI data using the WLS esti-mator.

(36)

in the presence of a time varying noise. It is clearly observed that the WLS method detects the false detection too. The Figure 3.16 shows the detection result of auditory cortical areas using the WLS estimator.

3.2.2.3 Bayesian estimator

From the Figures 3.13 and 3.15 the OLS and WLS estimators for the simu-lated data doesn’t give accurate result due to: i) The OLS estimator assumes i.i.d noise and does not consider time varying noise. ii) The WLS estimator requires accurate estimate of the noise covariance matrix. It uses the resid-ual of a traditional method (e.g. OLS) to estimate. It is difficult to perform an accurate estimation of the covariance matrix.

In the time varying noise model, the noise covariance matrix is a diag-onal matrix. This model is spatially and temporally non-stationary. This property is investigated by the Bayesian method. In the previous section, the Bayesian method was used to estimate the parameters w under the assumption of constant variance. In this section, the Bayesian method is extended to variance of noise changing with time, which gives an estimated noise structure closer to the true noise [26].

Figure 3.17: Detection results of simulated fMRI data using the Bayesian estimator.

(37)

Figure 3.18: Detection results of activation of auditory fMRI dataset using Bayesian estimator.

By observing Figure 3.17, it can be seen that the Bayesian estimator detects more true activations and less false activations in the presence of a time varying noise. By comparing Figures 3.13, 3.14, 3.15, 3.16, 3.17 and 3.18 it can be stated that the results of the Bayesian estimator is more ac-curate than the OLS and WLS methods.

In this section, the Bayesian estimator does not consider the drift in fMRI data. In the following section, this method is extended to drift model.

(38)

3.2.3 Drift model

The aim of fMRI data analysis is to detect the activated human brain re-gions based on inference drawn from the estimated parameter in GLM.

Drift is a low frequency signal that slowly varies across the complete pe-riod of the signal. The causes of drift are head movement during acquisition, local changes of magnetic field and internal physiological process like car-diac and respiration movements. Drift can be removed with a preprocessing procedure before statistical analysis, e.g. highpass filtering or the median method. It can also be removed by the drift model [28].

The measured fMRI signal contains three types of parameters: BOLD response, noise and drift. The GLM of measured data is

y = wb + f + , (3.11)

or based on description in Section 3.2.1.2,traditional way of GLM

y = Φw + , (3.12)

where y is the measured time series, f is the drift and is the noise with dimensions T × 1. Here, the parameter w is the scalar which represents the contribution of BOLD signal is to the measure the fMRI time series. The design matrix Φ = [b, f ] with dimension T × 2 and estimated parameter w = [w, 1].

The fMRI data exhibits fractional noise also. In the drift model, first ap-ply the discrete wavelet transform (DWT) to the fMRI data. The resulting wavelet coefficients at lower than fine scale (J0) are zero [29], because the

drift does not vary significantly below J0. The J0 value can be estimated

by using the Confidence Interval Criterion (CIC) as model selection crite-rion. It is more efficient than the Akaike Information Criterion (AIC) and the Schwartz Information Criterion (SIC). The activation parameters are ac-curately estimated by using the Bayesian estimator from the modified GLM. Figure 3.20 shows that the detection results of a simulated fMRI data by using Bayesian estimator. It can be effectively estimated and removed drift from fMRI data and detect true activations and less false detections. Figure 3.19 shows the result of activated detection. It is clearly seen that the activation of auditory cortical areas of the brain are successfully detected.

(39)

Figure 3.19: Detection result of Drift model on the auditory activation task.

Figure 3.20: Detection result of Drift model on the simulated fMRI dataset. In this chapter, methods are discussed on the detection of the activated

(40)

Chapter 4

Estimation

4.1 Estimation of the hemodynamic response

(HDR)

The HDR can be used to study the human brain functions. It exhibits temporal properties of human brain activities. Generally, the HDR can be estimated through a signal averaging procedure by assuming that the ISI is large [30], so that consecutive stimuli do not overlap. Overlapping of consecutive stimuli leads to bias results of the averaging methods. This method works only under the linear assumption case, i.e. ISI is larger than 4-6 sec, otherwise nonlinear phenomena are observed. In this chapter, I have described linear adaptive modelling and nonlinear method based on RBF neural networks to estimate the HDR and detection of brain activation. More over in this chapter describes the NARX to model the dynamics of fMRI data.

4.1.1 Adaptive spatiotemporal modeling

The fMRI data contains spatial and temporal information. In the adaptive spatial temporal modeling, spatial information is also included for estimat-ing the temporal activation. Because activated brain regions extend to sev-eral voxels, surrounding voxels are more likely to be activated compared to other voxels. This spatial information can be used to improve the SNR and detection accuracy [31].

(41)

4.1. ESTIMATION OF THE HEMODYNAMIC RESPONSE (HDR)

Figure 4.1: Schematic diagram of the spatial smoothing and temporal mod-elling.

Figure 4.1 shows a sketch of the spatial smoothing and temporal mod-elling. The spatial adaptive filter is defined as

d(n) =

L

X

i=0

wiyi(n) = wTy(n), (4.1)

where y(n) = [y0(n), y1(n), . . . , yL(n)]T, y0(n) is the nth sample of given

voxel fMRI time series and y1(n), y2(n) . . . yL(n) are L number of

neigh-boring voxels time serieses, w = [w0(n), w1(n), . . . , wL(n)]T is the

corre-sponding weight vector. The smoothed signal d(n) is the desired signal for temporal modelling. For activated voxel time series, the resulted smoothed signal is approximated to the ideal BOLD signal r(n) to the experimental stimuli, i.e.

d(n) = wTy(n) = r(n) + u(n), (4.2) where r(n) is the ideal BOLD signal which is the result of stimulus func-tion and canonical HDR and u(n) is the white noise. The output spatial smoothed signal d(n) is used as the desired signal for temporal modeling. The temporal modelling is defined as

d(n) =

P

X

m=0

s(n − m)hm+ (n) = hTs(n) + (n), (4.3)

where s(n) = [s(n), s(n − 1). . . s(n − P )]T _{is a vector corresponding to the}

stimulus function with a delay of P and h = [h0, h1. . . hP]T are the filter

coefficients.

(42)

E{(d(n) − hT_s(n))2_{}. The cost function is calculated as}

J = E{(r(n) − d(n))2} + E{(d(n) − hT_s(n))2_}

= E{(r(n) − wTy(n))2} + E{(wT_{y(n) − h}T_s(n))2_}. _(4.4)

For the cost function J , the optimum filter coefficients are obtained at minimum value of J . The least mean square (LMS) algorithm is used for finding minimum mean square error by updating the following equations,

e1(n) = r(n) − ˆwT(n)y(n), (4.5) e2(n) = wˆT(n)y(n) − ˆhT(n)s(n), (4.6) ˆ w(n + 1) = w(n) + 2µˆ 1e1(n)y(n) − 2µ1e2(n)y(n), (4.7) ˆ h(n + 1) = h(n) + 2µˆ 2e2(n)s(n). (4.8)

Figure 4.2: Schematic diagram of the adaptive spatio-temporal model. Figure 4.2 shows the structure of the adaptive spatio-temporal model and the estimation errors (e1(n) and e2(n))are feedback to spatial filter to

adjust the estimation coefficients. The estimation error (e2(n)) is fed back to

(43)

converges. The estimated optimum spatial filter weights ˆw approximates the spatial smoothing and the optimum temporal filter weights ˆw approximates HDR. The adaptive spatio-temporal model is tested on the simulated BOLD signal.

Figure 4.3: Estimation of HDR on the simulated fMRI signal. Figure 4.3 shows the simulated BOLD signal, simulated noisy signal and estimation of HDR.

(44)

4.1.2 Neural Network

As mentioned before, the fMRI signal exhibits linear and nonlinear prop-erties. It has been observed that when the duration of ISI is less than 4 to 6 sec, the event related fMRI signals exhibit nonlinear properties. The nonlinear properties of the BOLD signal can be expressed by the balloon model, including variables such as, blood flow, blood volume and oxygen concentration. It is a physiologically derived model which was introduced by Buxton et al [32]. It is a proper model for understanding the nonlinear mechanism underlying the neural activity and BOLD effect. The balloon model provides an approximate estimation of the parameters and it is not easy to estimate the parameters.

On the other hand, non-physiological models were developed for better estimation when the signal is highly non-linear. The Volterra series model is a dynamical input output system developed to estimate the HDR using first and second order kernels [16]. Higher the kernel orders, better the estimates. The number of parameters increases exponentially with higher order kernels. The performance depends upon the selected order of Volterra series. The lower order kernels may not capture the dynamic properties of the system as well as higher order kernels. Estimate of Volterra kernels by using the least squares (LS) method is difficult when the kernel order is high. However, the RBF neural network method efficiently estimates the HDR as well as it avoids singularity problem in the LS method [33].

Radial Basis Function (RBF) neural network

Neural networks is a powerful method for nonlinear modeling [34]. As shown in Figure 4.4, the present and past inputs of stimulus are applied to the neural network system and the output signal is the desired BOLD signal. The measured BOLD signal can be determined for given input stimulus s(n) as

ˆ

y(n) = FN N(s(n)), (4.9)

the output ˆy(n) is the nonlinear function of input s(n) and the nonlinear function is denoted by FN N. s(n) is a (P + 1) × 1 vector containing the present and past stimulus data with delays from 1 to P , s(n) = [s(n), s(n − 1) . . . s(n − p)]T.

(45)

Figure 4.4: Schematic diagram of the RBF neural network.

RBF network has a simple structure and is easy to implement compared to the multilayer perceptron of neural network. The output ˆy(n) is a linear combination of M radial basis functions defined as

ˆ y(n) = M X i=1 hiG(s(n), ci), (4.10)

where hiis the weighting coefficients and ci is the center of the radial basis

function. Commonly used basis functions are G(s(n), ci) = exp − 1 σ2 i ks(n) − cik2 . (4.11)

The center of a RBF is determined randomly from the data, and variance is selected according to the center spread. The aim of a RBF network is to find the weights and sum of minimized least square error. Regularization parameter is required for the stability of the solution. The weight vector is estimated by

(46)

For estimating the weight vector (h) and regularization, Bayesian learn-ing approach is used and is iteratively updated uslearn-ing the followlearn-ing equa-tions[35], Σ = β2(GTG + λI)−1, (4.13) ˆ h = 1 β2ΣG T_y, _(4.14) γ = M − λ β2trace(Σ), (4.15) β2 = ky − G ˆwk 2 N − γ , (4.16) λ = γβ 2 ˆ wT_w_ˆ, (4.17)

where β2 _{is the noise variance. Finally the BOLD signal is computed as}

ˆ

yb= Gˆh. (4.18)

Figure 4.5: Simulated BOLD signal, noisy BOLD signal and reconstructed BOLD signal by using RBF neural network.

(47)

Figure 4.5 shows the original simulated BOLD signal, simulated noisy signal with Gaussian noise and reconstructed BOLD signal by using RBF neural network method.

The detection of activated brain regions the reconstructed BOLD signal can be used to measures the activation of each voxel. The detection of activation index (R) is defined as

R = kˆybk ky − ˆybk

, (4.19)

where ˆyb is the reconstructed BOLD signal and y is the measured fMRI

time series in a voxel.

Figure 4.6: Detection results of simulated fMRI data with the RBF neural network model.

(48)

Figure 4.7: Detection results of auditory fMRI data using the RBF neural network method.

Figure 4.6 and Figure 4.7 show the result of activated brain regions of real fMRI data and simulated fMRI data using the RBF neural network method. By observing the results of auditory activation data, it can be clearly seen that the RBF neural network method detects the auditory cortical areas in the brain.

(49)

4.1.3 NARX model

The nonlinear autoregressive with exogenous inputs (NARX) neural network is a dynamical model for mapping the input-output dynamics of the BOLD signal. It is a powerful nonlinear model with faster convergence and better generalization ability [36].

Figure 4.8: Block diagram of the NARX model.

See Figure 4.8. In the NARX model, the measured fMRI signal y(n) and the input signal s(n) are applied as input to the NARX neural network. The input signal s(n) is applied to the NARX neural networks via q − 1 delays and the fMRI signal is applied to the NARX neural networks via p delays. The output estimated signal can be determined as

ˆ

yb(n) = FN N(y(n), s(n)) = FN N(x(n)). (4.20)

The estimated signal is a nonlinear transformation of the input signals s(n) and y(n). The input signal vector x(n) = [y(n), s(n)]Thas the size (p+q)×1. The output BOLD signal is estimated by the NARX model.

This method was developed through the RBF neural network, explained in Section 4.1.2. The BOLD signal is reconstructed after estimation of the weights. After reconstructing the BOLD signal, the following test can be

(50)

where ˆyb is the reconstructed BOLD signal and y is the measured fMRI

time series in a voxel.

Figure 4.9: Detection results of auditory fMRI data using the NARX model. From Figure 4.9, auditory activated areas can be detected by using the NARX model.

(51)

Chapter 5

Brain connectivity

The brain can be described as an organized structural network. Brain con-nectivity refers to anatomical concon-nectivity, functional concon-nectivity and ef-fective connectivity [38]. Anatomical connectivity is physical or structural connections between neurons or anatomical brain regions. Functional con-nectivity is a statistical concept, which represents the statistical dependency among different brain structures. Effective connectivity is the causal inter-action between different brain structures [39].

In this chapter, in the first section the functional connectivity is investi-gated using the nonlinear cross correlation analysis. In the second section, effective connectivity is investigated using the Granger causality.

5.1 Nonlinear cross correlation

The nonlinear cross-correlation coefficient describes the dependency of x on y without any specific relation between them. Here, x and y are the two time serieses. If the value of x is considered as function of y, then the value of y gives x and it can be predicted [40]. The nonlinear correlation coefficient is defined as, h2_y/x = PN k=1y 2_{(k) −}PN k=1(y(k) − f (x)) 2 PN k=1y2(k) , (5.1)

where f (x) is the approximation of the nonlinear method. In this work, f (x) is approximated with an RBF neural network. As explained in section 4.1.2, the present and past inputs of the time series x(n) is used as input for the RBF neural network. The output signal y(n) is the nonlinear function of input x(n) vector. The output signal can be expressed as [33]

(52)

5.2. GRANGER CAUSALITY

where y(n) = [y(n), y(n − 1) . . . y(n − p)]T _{is the time series y(n) vector of}

dimension (p + 1) × 1 formed by the present and past inputs with delay p. The output ˆy(n) is taken to be nonlinear function of the basis functions, it can be defined as ˆ y(n) = M X i=1 hiG(x(n), ci) (5.3)

where hi are the weighting coefficients, ci are the centers of the radial basis

functions and M is the hidden units of the RBF network. The weights can be estimated through the Bayesian learning procedure. When the weights are estimated, the functional approximation function can be expressed as

ˆ

y = Gˆh. (5.4)

The range of h2_y/x is between 0 to 1. The value 0 means that the signal y is independent of x. The value 1 represents that the signal y can be fully determined by the signal x [40]. The above process is applied to all voxels fMRI time series to the selected time series.

5.2 Granger causality

According to Norbert Wiener,“If the prediction of one time series can be improved by inclusion of the knowledge about other time series, then sec-ond time series is said to have causal influence on the first time series” [41]. However, this idea lacks the practical implementation. Later Granger proved practically, “If the error variance of one time series reduced by inclusion of lagged observation of second time series, then second time series has causal influence on first time series” [42].

Granger causality method is based on multiple linear regression for deter-mine whether a signal useful for forecasting another signal or not. According to Granger causality a time series x (or y) Granger causes a time series y (or x) if past values of x helps to predict the present value of y with better accuracy compared to considering its past values alone [43].

Jointly, two time series x and y can be represented by a bivariate regres-sive model x(n) = P X j=1 a2jx(n − j) + P X j=1 b2jy(n − j) + ε1(n), (5.5) y(n) = P X j=1 c2jx(n − j) + P X j=1 d2jy(n − j) + ε2(n), (5.6)

(53)

5.2. GRANGER CAUSALITY

In Equation (5.5) and (5.6), P is the maximum number of lagged ob-servations included in the model and n is the number of samples in signal. Generally P < n. The variables a, b, c and d are the coefficients of the model and ε1 and ε2 are the residual (prediction) errors for the each time

series. If the variance of ε1 (or ε2) is reduced by the inclusion of y (or x)

terms in the first (or second) equation, then it is said that y (or x) G-causes x (or y).

The magnitude of this interaction can be measured by using the log ratio of the prediction error variances for the restricted (R) and unrestricted (U) models,

fy→x= ln

var(ε1R)

var(ε1U)

. (5.7)

The restricted model can be defined as

x(n) = P X j=1 a1jx(n − j) + ε1R(n), (5.8) y(n) = P X j=1 d2jy(n − j) + ε2R(t), (5.9)

and the unrestricted model can be defined as

x(n) = P X j=1 a2jx(n − j) + P X j=1 b2jy(n − j) + ε1U(n), (5.10) y(n) = P X j=1 c2jx(n − j) + P X j=1 d2jy(n − j) + ε2U(n). (5.11) Model order

The estimation of multivariate autoregressive (MVAR) models requires the number of time-lags (P ) to include, i.e., the model orders. When the lags are lesser or too many these can lead to a poor representation or prob-lems in model estimation. A principled means to specify the model order to minimize the criteria that balance the variance accounted against the number of coefficients to be estimated. They can be implemented on two criteria’s, which are Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) [43]. The BIC is more useful in application to the neural systems and it is defined as

(54)

5.3. RESULTS

5.3 Results

5.3.1 Nonlinear cross correlation

Simulated fMRI data

First the proposed method was tested on simulated data. A simulated block design 3-D fMRI data set was investigated in this study. A block design BOLD signal was added at particular positions of the image, and then white Gaussian noise was added to the simulated signals with SNR values from 0.49 to 8.29 dB. Figure 5.1 shows the functional connected re-gions (blue color) with selected time series (yellow color).

Figure 5.1: Simulated image with activated brain regions (left) and func-tional connectivity of brain structures (blue pixels) of selected time series (yellow pixel) (right).

Real fMRI data

The auditory experimental data was used for investigation of functional con-nectivity of the brain structures. First, the activated brain regions was found by using one of the methods for detection of activated brain region, such as OLS, WLS or the Bayesian estimator (described in section 3.2). Then a time series voxel of interest was selected to map the functional connectivity. Figure 5.2 shows the reconstructed BOLD signal of functionally connected and unconnected voxel time serieses by using RBF neural network method. Figure 5.3 shows the functional connectivity of brain structures (dark-blue color) of the selected time series voxel (yellow color).

(55)

5.3. RESULTS

Figure 5.2: Time courses of functionally connected voxel time series and functionally unconnected time series with selected time series for the real fMRI data.

(56)

connec-5.3. RESULTS

the selected time series (yellow color voxel) is functionally connected with blue color voxel time series for auditory areas activation. The left auditory cortical areas is functionally connected with right auditory cortical areas during auditory activation.

Figure 5.4: Curve of nonlinear cross correlation values for different SNR values.

Figure 5.4 shows the nonlinear cross correlation values of two simulated time series for different noises. One of the time series is added with white Gaussian noise with different SNR values from −10 dB to 30 dB. The value of nonlinear cross correlation value without noise is 1. The nonlinear cross correlation is linear when SNR is changed from −10 dB to −2 dB, and then it varies nonlinearly until it reaches a maximum value (1), which is value for nonlinear cross correlation when there is no noise added.

5.3.2 Granger causality

Four simulated time series x1-x4 signals were generated by using

autore-gressive (AR) modeling. The relations between these signals are that x2 is

caused by x1 (x1 causes x2) and x3is caused by x4 (x4 causes x3), see

Fig-ure 5.5. These signals are added at particular positions of images containing Gaussian noise.

(57)

5.3. RESULTS

The simulated time series signals were x1(n) = 0.95 √ 2x1(n − 1) − 0.9025x1(n − 2) + w1(n), (5.13) x2(n) = 0.5x1(n − 2) + w2(n), (5.14) x3(n) = −0.4x4(n − 3) + w3(n), (5.15) x4(n) = 0.35x4(n − 2) + w4(n). (5.16)

Figure 5.5: Simulated image with causal relations of brain structures. The proposed method maps the causal relation of brain areas for selected fMRI signal time series voxels. From the Figure 5.6, the first figure shows that selected (yellow color) voxel time series is caused by red color voxels time series (red color voxel time series causes to yellow color voxel time series). In the second figure, the yellow color voxel time series (x3) is caused by red color voxel time series (x4) (red color voxel time series causes to yellow color voxel (x4) time series). The third figure, yellow color pixel not caused by any time series. By observing the results the Granger causality method gives proper result according to simulation in presence of Gaussian noise. From the Figure 5.7, in the real dataset the red color voxel time series or regions causes to yellow color voxel time series or region.

(58)

5.3. RESULTS

Figure 5.6: Causal relation of brain images: i) x2 (yellow pixel) caused by x1 time series (red pixels), ii) x3 (yellow pixel) caused by x4 (red pixels) time series iii) No causal relation with selected time series at the yellow pixel.

Figure 5.7: Real data set image with causal relation of brain structures (red pixels) with selected time series (yellow pixel).

Figure 5.8 shows the values of Granger causality interaction of two sim-ulated time series with different noises. One of the time series is added with white Gaussian noise with different SNR values ranging from −10 dB to 30 dB. For the given set of equations, the estimated Granger causality is 0.65 (i.e, without noise). When SNR is changed from −10 dB to 25 dB the Granger causality is nonlinear, later it reaches the maximum value 0.65, which is value when there is no noise.

(59)

5.3. RESULTS

(60)

Chapter 6

Conclusion and Future

work

6.1 Conclusion

Functional magnetic resonance imaging (fMRI) is a powerful technique to study the brain functions. It is MRI based technique to measure brain ac-tivity by detecting the changes in blood oxygen level. The analysis of fMRI data are challenging due to the complexity of the data and the presence of noise.

In the first part of thesis work, the developed Matlab toolbox provides buttons for display and loading fMRI data for the analysis. The sparse Bayesian learning algorithm is developed for construction of the design matrix in the GLM. This Bayesian method provides better brain activa-tion results and avoids multiple comparison problems in the t-test method. Through the comparison of ROC curve, the Bayesian estimator is robust than the t-test method. When the Bayesian approach is extended the time varying noise, the results are better than for the OLS and WLS methods. In the drift model, the Bayesian method was extended with the drift model. It successfully removes the drift from the fMRI data.

In the second part of this study, the linear and nonlinear methods for detection of brain activation and estimation of the HDR were implemented. The linear method in spatio temporal modeling is better choice when the Inter stimulus interval is large. The nonlinear RBF method is more flexible and computationally efficient. The NARX method developed to capture the dynamics of fMRI data.

The third part of study was mainly focused on brain connectivity meth-ods. A functional connectivity method is used to map the functionally

(61)

6.2. FUTURE WORK

connected regions of the brain structures by using nonlinear cross correla-tion analysis. This method could successfully map the funccorrela-tional connected regions in the brain. Moreover, effective connectivity is achieved by Granger causality approach. It determines changes in brain activation causal to the changes in other brain regions. The proposed method maps the causal in-teractions between brain structures. All these methods are investigated on both simulated and real datasets. The proposed Matlab toolbox in the thesis allows the investigator to study the brain functions and helps to do further research.

6.2 Future work

The developed toolbox provides results in axial view, the future work can be providing the results in 3D view. The Bayesian methods for the detection of brain activation were based on assumption of Gaussian noise. Normally, in MR imaging the Rician distribution of noise also exist. In the future the Bayesian method will be implemented in the presence of Rician distribution of noise.

The Granger causality was only applied for the bi-variate variables and the multivariate Granger causality was not performed. The relationship between more than two variables analysis with the Multichannel analysis methods such as the partial coherence, directed partial coherence need to be implemented. In this method, the fMRI serieses are assumed to be linear and stationary. The future work can be implemented multivariate Granger causality and nonlinear Granger causality.

(62)

Bibliography

[1] M. S. Gazzaniga, R. B. Ivry, and G. R. Mangun, Cognitive neuro-science : the biology of the mind, 2nd ed. W.W.Norton & Company, 2002, ch. 1.

[2] Martin A. Lindquist, The Statistical Analysis of fMRI Data, Statisti-cal Science 2008, Vol. 23, No. 4, pp. 439-464.

[3] S.A.Huettel, A.W.Song, and G.McCarthy, Functional Magnetic Res-onance Imaging. Sinauer Associates, Inc, 2004, ch.7.

[4] Govind B. Chavhan, MD, DNB , Paul S. Babyn, MD , Bejoy Thomas, MD , Manohar M. Shroff, MD , E. Mark Haacke, PhD, Principles, Techniques, and Applications of T2* based MR Imaging and Its Spe-cial Applications, pp. 1434-1448, 2009.

[5] S.A.Huettel, A.W.Song, and G.McCarthy, Functional Magnetic Res-onance Imaging. Sinauer Associates, Inc, 2004, ch.7 pp. 208-214. [6] P. Jezzard, P.M. Matthews, and S. M. Smith, Eds., Functional MRI:

an introduction to methods, Oxford University Press, 2001, ch. 8. [7] S.A.Huettel, A.W.Song, and G.McCarthy, Functional Magnetic

Res-onance Imaging. Sinauer Associates, Inc, 2004, ch.9, pp. 302-325. [8] SPM, By members and collaborators of the Wellcome Trust Centre for

Neuroimaging, Single subject epoch (block) auditory fMRI activation data, (http://www.fil.ion.ucl.ac.uk/spm/data/auditory/).

[9] S.A.Huettel, A.W.Song, and G.McCarthy, Functional Magnetic Res-onance Imaging. Sinauer Associates, Inc, 2004, ch.6, pp. 266-277. [10] SPM, By members and collaborators of the Wellcome Trust Centre

for Neuroimaging, SPM8 mauanl, ch. 28, pp. 209-216.

[11] Victoria L. Morgan, David R. Pickens, Steven L. Hartmann, and Ronald R. Price, Comparison of Functional MRI Image Realignment Tools Using a Computer-Generated Phantom. Magnetic Resonance in Medicine, vol. 46, pp. 510-514, 2001.

(63)

BIBLIOGRAPHY

[12] J. Talairach and P. Tournoux, A Co-planer Stereotaxic Atlas of the Human Brain. New York: Thieme Medical Publishers, 1988.

[13] A.W. Toga, Brian warping. Academic Press, 1999.

[14] K. Fristion, O. Josephs, E. Zarahn, A. Homes, S. Rouquette, and J. Poline, To smooth or not to smooth Bias and efficiency in fMRI time series analysis, Neuroimage, vol. 12, pp. 196-208,2000.

[15] K. Worsley. C. Liao, J. Aston, V.Petre, G.Duncan, F.Morales, and A. Evans, “A general statistical analysis for fMRI data,” Neuroimage, vol. 15, pp. 1-15, 2002.

[16] K. Friston, P. Fletche, O. Josephs, A. Holmes, M. Rugg and R. Turner, “Event-related fMRI:Characterising differential responses,” Neuroim-age, vol. 7, pp. 30-40, 1998.

[17] R. Frackowiak, K. Friston, C. Frith, R.Dolan, C. Price, S. Seki, J. Ash-burner, and W. Penny, Eds., HumanBrainF unction, 2nd ed. Elsevier Academic Press, 2003, ch. 37.

[18] J. Brosch, T. Talavage, J. Ulmer, and J. Nyenhuis, “Simulation of human respiration in fMRI with a mechanical model,” IEEE Trans-actions on Biomedical Engineering, vol 49, pp. 700-707, 2002. [19] A. Smith, B. Lewis, U. Ruttinmann, F. Ye, T. Sinnwell, Y. Yang,

J. Duyn, and J. Frank, “Investigation of low frequency drift in fMRI signal,” N euroimage, vol. 9, pp. 526-533, 1999.

[20] Huaien Luo and S. Puthusserypady, “A sparse Bayesian method for determination of flexible matrix for fMRI data analysis,” IEEE Trans-actions on Circuits and Systems, Part I Fundamental Theory and Ap-plications, Vol.52, No.12, pp.2699-2706, December 2005.

[21] Martina Pavlicov’a, Noel Cressie, and Thomas J. Santner, “Testing for Activation in Data from FMRI Experiments,” Journal of Data Science 4, pp. 275-289, 2006.

[22] Martin M. Monti, “Statistical Analysis of fMRI Time-Series: A Criti-cal Evaluation of the GLM Approach” Princeton University, Dept. of Psychology, April 18, 2006.

[23] R. Frackowiak, K. Friston, C. Frith, R. Dolan, C. Price, S. Zeki, J. Ashburner, and W. Penny, Eds., Human Brain Function, 2nd ed. El-sevier Academic Press, 2003, ch. 40.

(64)

eval-BIBLIOGRAPHY

[25] C.Long, E. Brown, C. Triantafyllou, I. Aharon, L.L. Wald, and V. Solo, “Nonstationary noise estimation in functional MRI,” NeuroIm-age, vol. 28, pp. 890-903, 2005.

[26] Huaien Luo and S. Puthusserypady, “fMRI data analysis with nonsta-tionary noise models: A Bayesian approach,” IEEE Transactions on Biomedical Engineering, Vol. 54, pp. 1621-1630, 2007.

[27] J. Diedrichsen and R. Shadmehr, “Detecting and adjusting for arti-facts in fMRI time series data,” N euroImage, vol 27, pp. 624-634, 2005.

[28] Huaien Luo and S. Puthusserypady, “Analysis of fMRI Data with Drift: Modified General Linear Model and Bayesian Estimator,” IEEE Transactions on Biomedical Engineering, vol. 55, pp. 1504-1511, 2008. 2008.

[29] F. Meyer, “Wavelet-based estimation of a semipara-metric generalized linear model of fmri time-series,” IEEET ransactionsonM edicalImaging, vol.22, no.3 pp.315-322,2003.

[30] A. Dale and R. Buckner, “Selective averaging of rapidly presented individual trails using fMRI,” Human Brain Mapping, vol. 5, pp. 329-340, 1997.

[31] Huaien Luo and S. Puthusserypady, “Adaptive Spatiotemporal Mod-elling and Estimation of the Event-related fMRI Responses,” Signal Processing, Vol. 87, pp. 2810-2822, 2007.

[32] R. Buxton, E. Wong, and L. Frank, “Dynamics of blood flow and oxy-genation changes during brain activation: the balloon model.” Mag-netic Resonance in Medicine, vol. 39, no. 6, pp.855-864, 1998. [33] Huaien Luo and S. Puthusserypady, “Estimation of the hemodynamic

response of the fMRI data using RBF neural networks,” IEEE Trans-actions on Biomedical Engineering, Vol. 54, pp. 1371-1381, 2007. [34] S. Haykin, Neural Networks, A comprehensive foundation. Prentice

Hall, 1999.

[35] E. Rank, “Application of Bayesian trained RBF networks to nonlinear time-series modeling,” Signal Processing, vol. 83, pp. 1393-1410, 2003. [36] Huaien Luo and S. Puthusserypady, “Spatio-temporal modeling and analysis of fMRI data using NARX neural network,” International Journal of Neural Systems, Vol.16, No.2, pp.139-149, 2006.

(65)

BIBLIOGRAPHY

[37] Mikail Rubinov, Olaf Sporns, “Complex network measures of brain connectivity: Uses and interpretations,” NeuroImage, vol. 52, pp. 1059-1069, 2010.

[38] Olaf Sporns, Indiana University, Bloomington, In, Brain Connectivity - Scholarpedia, 2007.

[39] Karl J. Friston, “Functional and Effective Connectivity in Neuroimag-ing: A Synthesis,” Human Brain Mapping, vol. 2, pp. 56-78, 1994. [40] Ernesto Pereda, Rodrigo Quian Quiroga, Joydeep Bhattacharya,

“Nonlinear multivariate analysis of neurophysiological signals,” Pro-gressive in Neurobiology. 77 pp. 1-7, 2005.

[41] Steven L. Bressler, Anil K. seth, “Wiener-Granger Causality: A well established methodology,” NeuroImage, vol. 58, pp. 323-329, 2011. [42] Andrea Brovelli, “Statistical Analysis of Single-Trail Granger

Causal-ity Spectra,” Computational and Mathematical Methods in Medicine, vol. 2012, pp. 1-10, 2012.

[43] Anil K.Seth, “A MATLAB toolbox for Granger Causal connectivity analysis.” Journal of Neuroscience Methods vol. 186 pp. 262-273, 2010.

(66)

(67)

“Kiran” — 2012/9/6 — 12:35 — page 57 — #67 Avdelning, Institution Division, Department Datum Date Spr˚ak Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats ¨Ovrig rapport

URL f¨or elektronisk version

ISBN ISRN

Serietitel och serienummer Title of series, numbering

ISSN

Link¨oping Studies in Science and Technology Thesis No. -12/4600--SE

Titel Title F¨orfattare Author Sammanfattning Abstract

Functional Magnetic Resonance Imaging (fMRI) is one of the best tech-niques for neuroimaging and have revolutionized the way to understand the brain functions. It measures the changes in the blood oxygen level-dependent (BOLD) signal which is related to the neuronal activity. Complexity of the data, presence of different types of noises and the massive amount of data makes the fMRI data analysis a challenging one. It demands efficient signal processing and statistical analysis methods. The inference of the analysis are used by the physicians, neurologists and researchers for better understanding of the brain functions.

The purpose of this study is to design a toolbox for fMRI data anal-ysis. It includes methods to detect the brain activity maps, estimation of the hemodynamic response (HDR) and the connectivity of the brain structures. This toolbox provides methods for detection of activated brain regions mea-sured with Bayesian estimator. Results are compared with the conventional methods such as t-test, ordinary least squares (OLS) and weighted least squares (WLS). Brain activation and HDR are estimated with linear adaptive model and nonlinear method based on radial basis function (RBF) neural network. Nonlinear autoregressive with exogenous inputs (NARX) neural network is developed to model the dynamics of the fMRI data. This toolbox also provides methods to brain connectivity such as functional connectivity

DIVISIONSHORT,

Dept. of Computer and Information Science 581 83 Link¨oping

05-09-2012

-LiU-Tek-Lic–2012:12/4600--SE

-A Matlab Toolbox for fMRI Data -Analysis: Detection, Estimation and Brain Connectivity

Kiran Kumar Budde

× ×

A Matlab Toolbox for fMRI Data Analysis: Detection, Estimation and Brain Connectivity

Master thesis

A Matlab Toolbox for fMRI Data

Analysis: Detection, Estimation and

Brain Connectivity

Kiran Kumar Budde

LiTH-ISY-EX--12/4600--SE

05-09-2012

Master thesis

A Matlab Toolbox for fMRI Data

Analysis: Detection, Estimation and

Brain Connectivity

Kiran Kumar Budde

LiTH-ISY-EX--12/4600--SE

05-09-2012

Abstract

Acknowledgement

Abbreviations

Contents

Chapter 1

Introduction

1.1

Introduction

1.2

Outline of the thesis

Chapter 2

Background

2.1

Functional magnetic resonance imaging

2.1.1

BOLD signal generation

2.1.2

Hemodynamic response (HDR)

2.1.3

fMRI Experimental design

2.1.4

fMRI experimental data used in this thesis

2.2

fMRI data analysis

2.2.1

Preprocessing

2.3

fMRI data modeling

2.3.1

Temporal modeling

2.3.2

BOLD Model

2.3.3

Noise and Drift

Chapter 3

Detection

3.1

Loading and visualization

3.1.1

Display

3.1.2

Loading real fMRI dataset

3.1.3

Loading simulated fMRI dataset

3.1.4

Loading stimulus

3.2

Detection of activated brain regions

3.2.1

Flexible design matrix

3.2.2

Nonstationary noise models

3.2.3

Drift model

Chapter 4

Estimation

4.1

Estimation of the hemodynamic response

(HDR)

4.1.1

Adaptive spatiotemporal modeling

4.1.2

Neural Network

4.1.3

NARX model