DECONVOLUTION OF NEURONAL SIGNAL FROM HEMODYNAMIC RESPONSE Martin Havlicek 1,2,4 ,Jiri Jan 1 , Milan Brazdil 3 , Vince D. Calhoun 2,4 1 Department of Biomedical Engineering, Brno University of Technology, Brno, Czech Republic 2 The Mind Research Network, Albuquerque, NM, USA 3 First Department of Neurology, St. Anne's University Hospital, Brno, Czech Republic 4 Dept of ECE, The University of New Mexico, Albuquerque, NM, USA havlicekmartin@gmail.com ABSTRACT In this paper we describe a deconvolution technique for obtaining an approximation of the neuronal signal from an observed hemodynamic response in fMRI data. Our approach, based on the Rauch-Tung-Striebel smoother for square-root cubature Kalman filter, enables us to accurately infer the hidden states, parameters, and the input of the dynamic system. Using a series of simulations we show in this paper that we are able to move beyond the limitation of a poorly sampled observation signal and estimate the true structure of underlying neuronal signal with significantly improved temporal resolution. Index Terms— cubature. Kalman, fMRI, deconvolution, smoother 1. INTRODUCTION In functional magnetic resonance imaging (fMRI) the observed hemodynamic response is an indirect measure of neuronal activation, i.e. the response is represented by changes in blood flow and blood oxygenation that follow after neuronal activation. Moreover, this complex relationship can be both subject and brain region specific, which makes identification of the true effective connectivity (directional influence) between different brain regions difficult. Fortunately, the hemodynamic model describing this relation has been characterized [1,2], which, under certain assumptions, allows an inversion of the hidden dynamic process to obtain an estimate of the neuronal activation given some observations. The Bayesian framework is the most commonly used method for the study of dynamic systems that allows inference on the hidden states, model parameters and the input. However, since the hemodynamic balloon model possesses strong nonlinear character, we need to utilize a nonlinear estimation procedure. The selection of suitable sub-optimal solutions to the recursive Bayesian estimation problem involved a trade-off between a) global optimality; and b) computational tractability and robust behavior. Taking this into account, we propose a deconvolution technique based on the recently introduced cubature Kalman filter (CKF) [3], which is the closest known direct approximation to the Bayesian filter and has been shown to outperform all other nonlinear filters in a Gaussian environment. Specifically, we apply an initial CKF step that is refined by a subsequent cubature Rauch-Tung-Striebel (RTS) smoother to simultaneously infer the hidden states, the input, and the parameters. 2. CONTINUOUS-DISCRETE DYNAMIC MODELS Nonlinear filtering problems are typically described by a state- space model consisting of a process equation and a measurement equation. In our case, the process equation is derived from the underlying physics of a continuous dynamic system, and expressed in the form of a set of differential equations. Whereas the measurements ݕare acquired in discrete times ሺ ݐൌ ͳǡʹǡ ǥ ǡ ሻ. Hence, we have a model with a continuous process equation and a discrete measurement equation in stochastic form:  ݔ ௧ ൌ ܐሺ ݔ ௧ ǡ ߠ ௧ ǡ ݑ ௧ ǡ ݐሻ ݐ൅ ܔሺ ݔ ௧ ǡ ݐሻ ߚ ௧ ݕ ௧ ൌ ሺ ݔ ௧ ǡ ߠ ௧ ǡ ݐሻ൅ ݎ ௧ ǡ (1) where ߠ ௧ represents unknown parameters of the equation of motion ܐand the measurement function , respectively; ݑ ௧ is the exogenous input that drives the hidden states; ݎ ௧ is Gaussian measurement noise; ܔሺ ݔ ௧ ǡ ݐሻ is a function of the states and time; and ߚ ௧ denotes Brownian motion with increments  ߚ ௧ being independent of the state process ݔ ௧ . The continuous time formulation of the stochastic differential equations (SDE) can be converted into a discrete-time analogue using e.g. local linearization (LL) scheme [4]: ݔ ௧ ൌ ܎ሺ ݔ ௧ଵ ǡ ߠ ௧ଵ ǡ ݑ ௧ଵ ሻ൅ ݍ ௧ଵ ǡ ݕ ௧ ൌ ሺ ݔ ௧ ǡ ߠ ௧ ሻ൅ ݎ ௧ ǡ (2) where ݍ ௧ is a Gaussian state noise vector of variance ܔ ଶ οݐ. In that case, the function ܎ሺǤ ሻ is evaluated through: ܎ሺǤ ሻ ൎ ݔ ௧ଵ ൅ ܬ ௫ ೟ ଵ ൫ ܬ ௫ ೟ οݐ൯െ ܫ൧ܐሺ ݔ ௧ଵ ǡ ߠ ௧ଵ ǡ ݑ ௧ଵ ሻǤ (3) Here ܬ ௫ ೟ is a Jacobian of ܐand ο ݐis the time interval between samples. The LL-scheme has demonstrated to improve the convergence and stability properties of conventional numerical integrators [4]. See [5] for its simple algebraic representation. 3. NONLINEAR MODEL IDENTIFICATION Parameter estimation, also referred to as a system identification can be regarded as a special case of general state estimation in which the parameters are treated as specific states. Parameter estimation involves determining nonlinear mapping: ݕ ௧ ൌ ሺ ݔ ௧ Ǣ ߠ ௧ ሻ, (4) where the nonlinear map ሺǤ ሻ is in our case the dynamic model ܎ሺǤ ሻ parameterized by the vector ߠ ௧ . The parameters ߠ ௧ correspond to a stationary process with identity state transition matrix, driven by an “artificial” process noise ݓ ௧ . Similarly, the input or the cause of the hidden states ݑ ௧ can be estimated as well. This is possible 617 978-1-4577-0539-7/11/$26.00 ©2011 IEEE ICASSP 2011