Nonlinear Bayesian Estimation of fMRI BOLD Signal under Non-Gaussian - - PowerPoint PPT Presentation

Nonlinear Bayesian Estimation of fMRI BOLD Signal under Non-Gaussian Noise

Ali Fahim Khan

Supervised by: Dr. Muhammad Shahzad Younis

• Dr. Kashif M. Rajpoot
• Dr. Khawar Khurshid
• Dr. Amir Ali Khan

Outline

 Introduction  Literature Review  Methodology  Experiments and Results  Conclusion and Discussion  Bibliography  Q&A Session  Appendix

Introduction

 fMRI-Functional Magnetic Resonance Imaging

 A window to the brain!  A non-invasive tool to study the neural activity

 Example: A blind individual reading Braille

Image courtesy http://cortivis.umh.es/overview.htm

Occipital Lobe Somatosensory Lobe

Applications of fMRI Making brain atlas (HCP) Brain disease diagnosis such as Alzheimer's (Sterling R., 2011). Training patients having brain illness to cure them (Birbaumer, et al., 2007) Lie detector (Langleben et al., 2005)

Introduction

• fMRI data is acquired using an fMRI

scanner

• Time series data generated from

voxels

• By studying this noisy data, one can

make inferences on brain activity

Photo Cour. devendradesmukh.blogspot.com

1 mm x 1 mm x 1.5 mm 7 mm x 7 mm x 10 mm

S.M Smith, “Overview of fMRI analysis”, The British Journal of Radiology

• It is this brain activity we

are interested in!

Literature Review

fMRI is most commonly performed using blood

• xygenation level-dependent

(BOLD) contrast (Ogawa et al.,1992)

(Douglas, 2001)

Methods to study BOLD signal

• Statistical Parametric Mapping (SPM) (Friston, 1995)
• Methods based on the Hemodynamic Model
• Hybrid methods such as Genetic Algorithms and simulated

annealing (Vakorin et al., 2007).

Literature Review

 The Hemodynamic Approach

1 1 1 ( 1)

1 1 ( ) ( 1) 1 ( ) 1 (1 ) 1

s f f

s u t s f f s v f v E q f v q E

 

    



                        

1 2 3 1 2 3

(1 ) (1 ) (1 ) 7 , 2, 2 0.2 q y V k q k k v v k E k k E                

First proposed by Buxton et al. 1998 Modified by Mandeville et al. 1999 Completed by Friston et al. 2000 Process Model Measurement Model

Literature Review

 Friston et al. first solved the hemodynamic model using Volterra

Kernels series (Friston et al. 2000)

 Later Gitelman et al. introduced Dynamic Causal Modeling that

linked different regions of the brain together (Gitelman et al. 2003). Does not include physiological noise

 Riera

et al. first included physiological noise to the hemodynamic model and performed blind deconvolution using Radial Basis Functions (Riera et al. 2004)

 Attempts to inverts the model using particle filters. However

Computationally demanding. (Johnston et al. 2008, Murray et al. 2008, Michah C. Chambers 2010)

Literature Review

 Hu et al. inverted the model using the SR-UKF (Hu et al. 2009)  Martin utilized the SR-CKF to perform blind deconvolution of

the Hemodynamic model (Martin, 2010)

 All of the methods above assumed both the process noise as

well as measurement noise to be Gaussian.

 Studies indicate Non-Gaussian noise in fMRI data

 Gamma distribution (Stephen et al., 2001)  Rician Distribution (Arnold et al., 2005)  Impulsive noise (Josephs at al., 2007)

Literature Review

 Gaussian Sum Filtering can tackle non-Gaussian noise. (Alspach

and Sorenson 1972). However computationally expensive

 To alleviate the computational problem, (Plataniotis et al. 1997)

proposed the Adaptive-GSF (AGSF). He and colleagues demonstrated its use in narrowband inference in presence of non-Gaussian noise (Plataniotis et al. 2000)

 (Miroslav et al. ,2005) proposed the Sigma point GSF (SPGSF)

having bank of SR-UKF and applied it to the non-Gaussian noise distribution problem.

Problem Statement

The problem boils down to finding the estimates of given a noisy times series data .

k

x

k

y

Methodology

 Modified Adaptive Gaussian Sum Filter

EKF working in parallel each tuned to a component of the noise

• f interest. The posterior probability distribution is subsequently

collapsed to yield single Gaussian term

• 0.02

0.02 0.04 0.06 0.08 0.1 0.12 5 10 15 20 25 30 35 Amplitude Difference Probability Density Gaussian pdf 1 Gaussian pdf 2 Gaussian pdf 3 Gaussian pdf 2 Composite non-Gaussian pdf

Methodology

 Sigma Point Gaussian Sum Filter (SPGSF)

Quite similar to GSF. But here multiple SR-UKF work in parallel each tuned to a specific Gaussian term constituting the Gaussian mixture needed to obtain the closed loop Bayesian Recursive Relation

Experiments

50 100 150

• 0.5

0.5 1

Time(Seconds) s

Input Stimulus 50 100 150 1 2 3

f Time(Seconds)

50 100 150 0.5 1 1.5

v Time(Seconds)

50 100 150 0.5 1 1.5

q Time(Seconds)

50 100 150

• 0.05

0.05

Clean BOLD Time(Seconds) Voxel Data Synthesis

50 100 150

• 0.1

0.1

Noisy BOLD Time(Seconds)

Table from (Friston et al. 2000)

TR=1.2s 13 seconds ON=1 and 13 OFF=0 periodic stimulus for 150 Seconds [0 1 1 1] Initial states Block design Exp. parameters

MAGSF vs EKF (Synthetic Data)

 Impulsive noise (e-mixture)

where is the mixing parameter and varies from 0.01 to 0.1 and the ratio of is usually in the range of 10 to 10,000 (vatola ,1984)

 Due to flexible nature, e-mixture can be used to model several

non-Gaussian distributions

1 1 2 2

2 2

~ (1 ) ( , ) ( , )

n n

w e N eN      

(0,1)

e

2 1

/  

1 1 2 2 1 2 1 2

2 2 2 2

~ (1 ) ( , ) ( , ) 0.01,

0.01, 0.02, 0.0001, 0.05

n n n n

w e N eN e           

   

0.0223 (9) eye 

v

Q .

Process Noise Measurement Noise

2 1

ˆ ( )

k true i

RMSE k



   x x

Performance Criteria

SPGSF vs SR-UKF (Synthetic Data)

 Impulsive noise (e-mixture), same as used for MAGSF  Gamma noise

0.0223 (9) eye 

v

Q .

Process Noise Measurement Noise

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 5 10 15 20 25 30 35 X Frequency

( ) 0.048 (0.00299,3.07 7) 0.318 (0.018,1.93 5) 0.109 (0.334,5.31 5) 0.525 (0.011,1.06 5)

K

p w N e N e N e N e            

20 40 60 80 100 120 140 0.01 0.02 0.03 0.04 0.05

Sample Number Error Amplitude

Error Histogram Error Realization

1 2 1 2

2 2

0.01,

0.02, 0.03, 0.0001, 0.05

n n

e     

   

Error Amplitude Frequency

State Estimation Results (Synthetic Data)

MAGSF vs EKF under Impulsive Noise SPGSF vs SRUKF under Impulsive Noise SPGSF vs SRUKF under Gamma Noise

Parameter Estimation Results(Synthetic Data)

MAGSF vs EKF under Impulsive Noise SPGSF vs SR-UKF under Impulsive Noise SPGSF vs SR-UKF under Gamma Noise

RMSE Evolution with Time(Synt. Data)

MAGSF vs EKF under Impulsive Noise SPGSF vs SR-UKF under Impulsive Noise SPGSF vs SR-UKF under Gamma Noise

2 1

ˆ ( )

k true i

RMSE k



   x x

Performance Criteria

Experiment using Real fMRI Data

 Preprocessing was

performed using SPM8 and brain activation map was generated.

 Voxels/ROI with highest

activation were selected and their time series extracted

 Marsbar Used for

Voxel/ROI management Block design experiment . Right hand clenched then left hand periodically TR=1.92s 13 seconds ON=1 and 13 OFF=0 Data taken from mccauslandcenter.sc.edu

MAGSF vs EKF, Impulsive Noise

State Estimates Parameter Estimates

SPGSF vs SRUKF, Impulsive Noise

State Estimates Parameter Estimates

Summary of Results

Results of our study using real data Results of previously reported studies

Conclusion and Discussion

 A filter designed to work under Gaussian noise may perform

poor if made to work under non-Gaussian noise.

 We proposed a novel filter (MAGSF) and applied it together

with SPGSF to the framework of the hemodynamic model. Our results match with those presented in earlier studies.

 Both are global filters and are less susceptible to getting stuck

at a local minima

 The parallel architecture enables noise of arbitrary distribution

to be handled by the filters.

 The proposed filter can be applied to a number of applications

such as sensor fusion, radar tracking applications, etc.

 Since the filters are recursive, they can be implemented on

hardware with ease.

Contributions to the Scientific Community

Proposed a novel optimal filter and applied it to the problem of hemodynamic model under non-Gaussian noise environment “Nonlinear Bayesian Estimation of BOLD Signal Under Non-Gaussian Noise” (Accepted and in Press) by Computational and Mathematical Methods in Medicine. Impact Factor 1.018 Applied SPGSF to the problem of hemodynamic model under non-Gaussian noise “Nonlinear Bayesian Estimation of BOLD Signal Under Non-Gaussian Noise Using Sigma Point Gaussian Sum Filter” (Under Review) by EURASIP Journal of Advances in Signal Processing. Impact Factor 0.81

Bibliography



Human Connectome Project, http://www.humanconnectomeproject.org



Sterling R., “Potential of functional MRI as a biomarker in early Alzheimer's disease,” Neurobiol Aging. 2011 Dec;32 Suppl 1:S37-43



Birbaumer N, Cohen LG., “Brain-computer interfaces: communication and restoration of movement in paralysis,” J Physiol. 2007 Mar 15;579(Pt 3):621-36



Langleben DD, Loughead JW, Bilker WB, Ruparel K, Childress AR, Busch SI, Gur RC. Telling truth from lie in individual subjects with fast event-related fMRI, Hum Brain Mapp. 2005 Dec;26(4):262-72.



Douglas, “A primer on MRI and FMRI”, PhD thesis, 2001



OGAWA,S.,TANK,D.,MENON,R.,ELLERMAN,J.,KIM,S., MERKLE,H.andUGURBIL, K. (1992). Intrinsic signal changes accompanying sensory simulation: Functional brain mapping and magnetic resonance

• imaging. Proc. Nat. Acad. Sci. 895951–5955.



• K. J. Friston, A.P

. Holmes, K.J. Worsley, J-B. Poline, C.D. Frith, and R.S.J. FrackoWiak. “Statistical parametric maps in functional imaging: A general linear approach,” Human Brain Mapping, vol. 2, no. 4,

• pp. 189-210, 1995.



Friston, K.J., Mechelli, A., Turner, R., Price, C.J., 2000. Nonlinear responses in fMRI: the balloon model, Volterra kernels, and other hemodynamics. NeuroImage 12, 466–477.



Vasily A. Vakorin, Olga O. Krakovska, Ron Borowsky, and Gordon E. Sarty. Inferring neural activity from BOLD signals through nonlinear optimization. NeuroImage, 38(2):248–60, November 2007



Buxton, R.B., Wong, E.C., Frank, L.R., “Dynamics of blood flow and oxygenation changes during brain activation: the balloon model,” Magnetic Resonance Med. 39, pp. 855-864, 1998.



Mandeville, J.B., Marota, J.J.A., Ayata, C., Zaharchuk, G., Moskowitz, M.A., Rosen, B.R., Weisskoff, R.M. “Evidence of cerebrovascular postarteriole windkessel with delayed compliance,” Journal of Cerebral Blood Flow Metabolism. 19 (6), pp. 679–689, 1999.

Bibliography



K.S. Vatola, “Threshold detection in narrowband non-Gaussian noise,” IEEE Transaction on Communications, vol. COM-32, pp. 134-139, 1984.



Gitelman, D.R., Penny, W.D., Ashburner, J., Friston, K.J., “Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution,” Neuroimage 19, pp. 200-207, 2003.



Riera, J., Watanabe, J., Kazuki, I., Naoki, M., Aubert, E., Ozaki, T., Kawashima, R., “A state- space model of the hemodynamic approach: nonlinear filtering of BOLD signals,” NeuroImage 21, pp. 547-567, 2004.



Johnston, L.A., Duff, E., Mareels, I., Egan, G.F., “Nonlinear estimation of the BOLD signal,” NeuroImage 40, pp. 504-514, 2008.



Murray, L., Storkey, A., “Continuous time particle filtering for fMRI,” Advances in Neural Information Processing Systems, 20, pp. 1049–1068, 2008.



Michah C. Chambers, “Full Brain Blood-Oxygen-Level-Dependent Signal Parameter Estimation Using Particle Filters”, MS-Thesis, Virginia Polytechnic Institute and State University, 2010.



Zhenghui Hu, Xiaohu Zhao, Huafeng Liu, and Pengcheng Shi. “Nonlinear Analysis of the BOLD Signal,” EURASIP Journal on Advances in Signal Processing, vol. 2009:Article ID 215409, 2009.

Bibliography

 M. Havlicek, K. Friston, J. Jan, M. Brazdil, and

V.D. Calhoun, "Dynamic Modeling of Neuronal Responses in fMRI using Cubature Kalman Filtering, " NeuroImage, vol. 56,

• pp. 2109-2128, 2011.

 Stephen J. Hanson, Benjamin M. Bly, “The Distribution of BOLD Susceptibility effects

in the Brain is Non-Gaussian,” Neuroreport.;12(9):1971-7, 2001

 Arnold J. D. Dekker, Jan Sijbers, “Implications of the Rician distribution for fMRI

generalized likelihood ratio tests,” Magnetic Resonance Imaging, 23, pp. 953–959, 2005

 O. Josephs, N. Weiskopf, R. Deichmann, “Detection and correction of spikes in fMRI

data,” Proceedings International Society Magnetic Resonance Med., 15, 2007

 H. Sorenson, D.L. Alspach, “Recursive Bayesian estimation using Gaussain sum,”

Automatica, vol. 7, pp.465-479, 1971.

 K.N. Plataniotis, D. Androutsos, A.N.

Venetsanopoulos, “Nonlinear filtering of non- Gaussian noise,” Journal of Intelligent and Robotic Systems, vol. 19, pp.207-231, 1997.

 K.N. Plataniotis, A.N.

Venetsanopoulos, “State Estimation in the Presence of non- Gaussian Noise,” AS-SPCC, 2000.

Q&A Session

Thank you!

Appendix

Noise Analysis

Noise Analysis done on both the active and inactive voxels

50 100 150 200 250 300 0.2 0.4 0.6 0.8 1 Sample Number) Normalized Amplitude

Noise Analysis

Most Active Voxels

50 100 150 200 250 300 0.88 0.9 0.92 0.94 0.96 0.98 1 Sample Number) Normalized Amplitude

Noise Analysis

Least Active Voxels

50 100 150 200 250 300 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Sample Number) Normalized Amplitude

Noise Analysis

Correlation Matrix between normalized adjacent most activated voxel time series.

V1 to V5 represents voxel 1 to voxel 5 respectively.

Voxel V1 V2 V3 V4 V5 V1 1.0000 0.8839 0.5406 0.8977 0.3716 V2 0.8839 1.0000 0.4855 0.8424 0.4197 V3 0.5406 0.4855 1.0000 0.6149 0.7723 V4 0.8977 0.8424 0.6149 1.0000 0.5038 V5 0.3716 0.4197 0.7723 0.5038 1.000

Noise Analysis

Correlation Matrix between normalized adjacent most activated voxel and least activated voxel time series.

Voxel V1 V2 V3 V4 V5 V6 V7 V1 1.0000 0.7135 0.0873 0.0007

• 0.0356
• 0.2654

0.1059 V2 0.7135 1.0000 0.3408 0.5555

• 0.0236
• 0.2855

0.0331 V3 0.0873 0.3408 1.0000 0.2797 0.0736

• 0.0825

0.0579 V4 0.0007 0.5555 0.2797 1.0000 0.4637 0.0946 0.1714 V5

• 0.0356
• 0.0236

0.0736 0.4637 1.0000 0.1036 0.5672 V6

• 0.2654
• 0.2855
• 0.0825

0.0946 0.1036 1.0000 0.1367 V7 0.1059 0.0331 0.0579 0.1714 0.5672 0.1367 1.0000

V1 to V7 represents voxel 1 to voxel 7 respectively.

Noise Analysis

Correlation Matrix between normalized adjacent least activated Voxel time series.

V1 to V6 represents voxel 1 to voxel 6 respectively.

Voxel V1 V2 V3 V4 V5 V6 V1 1.0000

• 0.0286

0.0744 0.1054

• 0.0058
• 0.0499

V2

• 0.0286

1.0000 0.7135 0.0007

• 0.0356

0.1059 V3 0.0744 0.7135 1.0000 0.5555

• 0.0236

0.0331 V4 0.1054 0.0007 0.5555 1.0000 0.4637 0.1714 V5

• 0.0058
• 0.0356
• 0.0236

0.4637 1.0000 0.5672 V6

• 0.0499

0.1059 0.0331 0.1714 0.5672 1.0000

Noise Analysis

Difference between two time series of voxels yields noise

50 100 150 200 250 300 350

• 0.1
• 0.05

0.05 0.1 0.15 Sample Number Signal Amplitude Voxel Adjacent Voxel 50 100 150 200 250 300 350 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Sample Number) Amplitude Difference 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 5 10 15 20 25 30 Amplitude Difference Density data Normal Distribution Gamma Distribution

• 3
• 2
• 1

1 2 3

• 0.01

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

Standard Normal Quartiles Amplitude Difference

Noise Analysis

Q-Q plot to test Gaussianity

• 3
• 2
• 1

• 0.1
• 0.08
• 0.06
• 0.04
• 0.02

Standard Normal Quartiles Amplitude Difference

• 3
• 2
• 1

• 0.2
• 0.15
• 0.1
• 0.05

Standard Normal Quartiles Amplitude Difference

• 3
• 2
• 1

• 0.06
• 0.04
• 0.02

Standard Normal Quartiles Amplitude Difference

• 3
• 2
• 1

• 0.12
• 0.1
• 0.08
• 0.06
• 0.04
• 0.02

Standard Normal Quartiles Amplitude Difference

• 3
• 2
• 1

• 0.06
• 0.05
• 0.04
• 0.03
• 0.02
• 0.01

Standard Normal Quartiles Amplitude Difference

• 3
• 2
• 1

• 0.06
• 0.04
• 0.02

Standard Normal Quartiles Amplitude Difference