Computational and Mathematical Methods in Medicine

Volume 2015, Article ID 389875, 10 pages

http://dx.doi.org/10.1155/2015/389875

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

^{1}School of Electrical Engineering and Computer Science, National University of Sciences and Technology, Islamabad 44000, Pakistan^{2}Department of Electrical Engineering, Institute of Space Technology, Islamabad 44000, Pakistan^{3}Department of Telecommunication Engineering, University of Engineering and Technology, Taxila 47080, Pakistan

Received 2 August 2014; Revised 13 October 2014; Accepted 20 October 2014

Academic Editor: Giancarlo Ferrigno

Copyright © 2015 Ali Fahim Khan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Modeling the blood oxygenation level dependent (BOLD) signal has been a subject of study for over a decade in the neuroimaging community. Inspired from fluid dynamics, the hemodynamic model provides a plausible yet convincing interpretation of the BOLD signal by amalgamating effects of dynamic physiological changes in blood oxygenation, cerebral blood flow and volume. The nonautonomous, nonlinear set of differential equations of the hemodynamic model constitutes the process model while the weighted nonlinear sum of the physiological variables forms the measurement model. Plagued by various noise sources, the time series fMRI measurement data is mostly assumed to be affected by additive Gaussian noise. Though more feasible, the assumption may cause the designed filter to perform poorly if made to work under non-Gaussian environment. In this paper, we present a data assimilation scheme that assumes additive non-Gaussian noise, namely, the -mixture noise, affecting the measurements. The proposed filter MAGSF and the celebrated EKF are put to test by performing joint optimal Bayesian filtering to estimate both the states and parameters governing the hemodynamic model under non-Gaussian environment. Analyses using both the synthetic and real data reveal superior performance of the MAGSF as compared to EKF.

#### 1. Introduction

The understanding of neuronal activation of brain, associated with a particular cognitive event, has caused great interest in the neuroimaging community. Functional magnetic resonance imaging (fMRI) modality offers an attractive mechanism to get an insight into the hemodynamic changes occurring inside the brain by exploiting the BOLD signal with high spatial resolution. The de facto method to analyze BOLD signal still remains the statistical parametric (SPM) technique developed by Friston et al. [1]. It attempts to establish a linear relationship between the observed BOLD signal and the underlying physiological processes inside the brain by convolving so-called Hemodynamic Response Functions (HRFs) with a specific design matrix corresponding to the experimental stimulus. Unfortunately, HRFs vary considerably from one part of the brain to the other and from person to person [2, 3]. An approach to battle the variation of HRF across the brain was presented by Gandolla et al. [4] where including the HRF spatial and temporal derivates were able to better capture the HRF inhomogeneities across the brain.

Although neurons in the brain exist at microscopic level, fMRI studies are made at macroscopic level by analyzing composite data generated from a plethora of neurons residing in small regions in the brain called voxels. Attempts have been made to model the BOLD signal by a set of nonlinear, nonautonomous differential equations linking hemodynamic changes with a set of psychological variables. Balloon or simply the hemodynamic model is one such example wherein extending the idea of distending venous chamber [5] and the Windkessel theory [6], Friston et al. establish the neurovascular coupling of physiological variables, namely, neuronal stimulus signal, cerebral blood volume (CBV), cerebral blood flow (CBF), and total deoxyhemoglobin content (dHb) [7]. In an fMRI experiment, each voxel produces a noisy time series data constituting the observation data which is used to invert the hemodynamic model to estimate the states as well as the parameters. Friston et al. calculated the hemodynamic model parameters using a Volterra Kernel series approximation to characterize the hemodynamic response [7]. Later they introduced a Bayesian framework free of Volterra Kernels [8]. Coupled brain regions were included in the model later on [9]. Such type of model inversion is known as dynamic causal modeling (DCM) and is extensively used to study brain connectivity. Though popular, the DCM does not incorporate physiological noise. This assumption is impoverished as studies contend the presence of significant endogenous fluctuations in neuronal activity [10]. Riera et al. went a step further and included physiological noise in their state space stochastic model. More interestingly, they performed blind deconvolution by using radial basis functions (RBS) to estimate the neuronal input signal [11]. Their implementation utilized the local linearization filter (LLF) [12] that closely resembled the celebrated extended Kalman filter [13] for continuous dynamic systems.

To accommodate for the high nonlinearities in the hemodynamic model, Johnston et al., Murray and Storkey and Michah employed the particle filter for model inversion [14–16]. Unlike the Kalman filter family which is a local filter, particle filter is a global filter. It handles the nonlinearities well and offers superior performance but at the cost of high computation. Utilizing the nonlinearity combating capability filter, Hu et al. utilized the square root unscented Kalman (SRUKF) filter to perform joint state and parameter estimation [17]. Later, they proposed a dual UKF filter wherein the states and parameters were estimated separately. The scheme eliminated artificially induced and unwanted covariance between the state and parameter estimates and also offered faster convergence [18]. Havlicek et al. utilized the square root cubature Kalman filter (SRCKF) to perform joint estimation of the physiological states and parameters of the balloon model [19]. A similar technique presented by Hu et al. [20] catered for low frequency drift often found in fMRI scans. Out of the above stated techniques, only Riera and Martin performed blind deconvolution; that is, they estimated not only the states and parameters but also the input signal. Another set of techniques other than Bayesian have been introduced by Friston named dynamic expectation maximization (DEM) and variational filtering [21].

The techniques mentioned above have contributed significantly to the understanding and quantification of the physiological processes responsible for neuronal activation. Unfortunately, these approaches mostly assumed both the process noise driving the dynamic model and the measurement noise affecting the measured data to be Gaussian in nature. Although the assumption offers the attraction of simplicity of implementation at face value, yet it is pernicious. Not surprisingly, real world systems are often plagued by non-Gaussian noises such as impulsive noise, Rayleigh noise, Levy noise, and Cauchy noise and so forth. In some studies, fMRI data has been shown to follow heavy tailed gamma distribution [22] and Rician distribution [23]. Impulsive noise is another common non-Gaussian noise source that is common in fMRI time series [24]. In such an environment, the filters designed for Gaussian noise are susceptible to poor performance and even divergence.

Several strategies may be employed to handle non-Gaussianity. An example is Middleton classes A, B, and C models [25, 26]. Another potent technique presented by Sorenson and Alspach, known as the Gaussian sum filter (GSF), considered both the state noise and the measurement noise to be non-Gaussian [27]. Their assumption was that a series of Gaussian distributions could be combined together to represent the noise sequences. They, however, used a finite truncated series of Gaussian terms to yield an optimal minimum mean square estimator (MMSE) filter. Thus, the GSF had a bank of Kalman filters working in parallel, each filter matched to a single term of the Gaussian sum. A major setback of this approach was the growing computational burden at each iteration step. Thus the method is not computationally feasible in several applications. To mitigate the computational cost, Masreliez and Martin proposed a filter based on the so-called score-function [28, 29]. Although less computationally expensive than the GSF, the proposed strategy required problem dependent nonlinearities and still a computationally expensive procedure to develop the score function.

Adaptive GSF is a filter that is an approximate version of the GSF and can be used without the problem of increasing computational burden [30]. The filter employs a Bayesian learning method to transform a mixture of non-Gaussian terms into a single Gaussian term. Plataniotis and Venetsanopoulos demonstrated the use of AGSF in the application of narrowband interference suppression under non-Gaussian noise [31].

In this paper we attempt to invert the Balloon model in the presence of non-Gaussian measurement noise while still keeping the process noise to be Gaussian. The noise chosen is the -mixture impulsive noise which is a form of Gaussian mixture. A modified adaptive Gaussian sum filter having EKF (MAGSF) working in parallel is employed to jointly estimate the states and parameters of the hemodynamic model. For comparison, the MAGSF is compared with EKF. Subsequent analyses show superior performance of MAGSF as compared to the EKF in non-Gaussian noise.

#### 2. Materials and Methods

##### 2.1. Hemodynamic Model

The hemodynamic model describes a mechanism that links the synaptic activity to the observed time series fMRI BOLD signal in a voxel. Although several enhancements to the original model have been proposed [32–34], the model initially proposed by Buxton et al., Mandeville et al. and later completed by Friston et al. is sufficient [35]. Figure 1 shows an example of the hemodynamic approach. A human test subject is subjected to a periodic visual stimulus. Consequently, the brain generates a flow inducing electric signal that causes blood to flow into a venous compartment fed by blood capillaries inside a region in the brain. This inflatable venous compartment can be further modeled as an input state output system having compartment volume and deoxyhemoglobin (dHb) content as the two state variables. The input to this system is blood flow and the resultant output is the BOLD signal which is the observed signal obtained by the fMRI acquisition equipment. As the blood flow increases, the venous “balloon” is inflated resulting in the expulsion of deoxygenated blood at a greater rate. However, this expulsion is not sufficient enough to counteract the influx of oxygenated blood to the venous compartment resulting in the “early dip” in the BOLD signal. A peak in the signal is observed with a time lag of about 4–6 seconds following the onset of neural activity. Once the distended balloon relaxes, the reduced compartment volume and blood dilution result in another dip following the peak known as the poststimulus undershoot [7].