Reliable and Efficient Approach of BOLD Signal with Dual Kalman Filtering

Abstract

By introducing the conflicting effects of dynamic changes in blood flow, volume, and blood oxygenation, Balloon model provides a biomechanical compelling interpretation of the BOLD signal. In order to obtain optimal estimates for both the states and parameters involved in this model, a joint filtering (estimate) method has been widely used. However, it is flawed in several aspects (i) Correlation or interaction between the states and parameters is incorporated despite its nonexistence in biophysical reality. (ii) A joint representation for states and parameters necessarily means the large dimension of state space and will in turn lead to huge numerical cost in implementation. Given this knowledge, a dual filtering approach is proposed and demonstrated in this paper as a highly competent alternative, which can not only provide more reliable estimates, but also in a more efficient way. The two approaches in our discussion will be based on unscented Kalman filter, which has become the algorithm of choice in numerous nonlinear estimation and machine learning applications.

1. Introduction

A thorough understanding of the dynamic relationship between cerebral blood flow (CBF), cerebral blood volume (CBV), and the blood oxygenation level dependent (BOLD) signal is essential for the physiological interpretation of fMRI activation data. The Balloon model described by Buxton et al. (1998) [1] is the first biomechanical plausible model to expound this relationship: increasing the flow (or perfusion rate) generally leads to dilution of venous deoxyhemoglobin (dHb), reducing the tendency of the blood to attenuate the magnetic resonance signal. The resultant increase in signal intensity is referred to as the BOLD response [2]. It is by extending this model to cover the dynamic coupling between CBF and synaptic activity, more sophisticated physiological realities are incorporated, for example, oxygen metabolism dynamics, both intra- and extravascular signal [3, 4], and more intricate models obtained.

The Balloon model is an input-state-output model with three state variables: blood flow, and volume, deoxyhemoglobin content and several biologically reasonable parameters. The problems of state estimation and parameter estimation (sometimes referred to as system identification or machine learning) associated with Balloon model are often formulated in a state-space representation, where Balloon model serves as a set of continuous-time system equations to describe the hemodynamic process. The equations are nonlinear, corresponding to the fact that Balloon model is one of the numerous nonlinear approaches to characterizing evoked hemodynamic response in fMRI.

Several work has utilized the Balloon model or its enhanced versions in the analysis of fMRI response. Some approaches, including expectation maximization (EM) [5, 6] and maximum likelihood [7], model the BOLD observation as deterministic hemodynamic process. However, a limitation of these methods is that they can only deal with measurement noise. Many promising approaches to the dual estimation problem belong to filtering algorithm that is able to account for both the physiological and measurement noise. Riera et al. [8] addressed the data assimilation problem in an extended Kalman filter (EKF) strategy. As EKF might lead to the problem of divergence due to linearized approximation, Johnston et al. proposed particle filter [9, 10] to avoid the flaw of linearization. Moreover, Hu et al. [11–13] employed unscented Kalman filter (UKF) that also outperforms EKF in terms of estimation error but with roughly the same computational cost. Most recently, Friston et al. described variational filtering to optimize the approximation of posterior density on hidden model variables, while accumulating sufficient statistics to optimize the conditional densities of parameters and precision [14, 15].

The approaches mentioned above have greatly improved our ability to explore, and above all, to quantify the physiological mechanism involved in neural activation. However, they still have palpable defects. It is noteworthy that many of them actually are in the spirit of joint filtering, in which the underlying states and parameters are concatenated into a single higher dimensional joint state space, a filter runs for estimating both the states and parameters. Despite its straightforwardness in theory and convenience in implementation, the weakness of joint filtering is obvious. The objective of this paper is to introduce and develop an estimator equally concise but with higher performance—dual filtering. We will demonstrate the advantage of dual filter from two aspects, by the example of dual UKF versus joint UKF. (1) In terms of Balloon model, there is no inherent biophysical correlation between the states and parameters. By treating them separately, dual filtering can avoid undesired transaction between them. (2) Larger dimension of state-space vector implies much more computational expense. Specifically, computational complexity for general state-space problems is 𝒪(𝐿³) [16]. Although the frequency of predict-update cycle required by dual filter is the twice of that required by joint filter, dual estimate is much more computational efficient.

2. Materials and Methods

2.1. Hemodynamic Model

Balloon model describes the coupled kinetic changes from synaptic activity to the fMRI BOLD signal at a given region. This model has been extended by Friston et al. (2000) [5] to include the effects of external inputs to an autoregulated vasodilatory signal, assuming that the relationship between evoked neural activity and blood flow is linear. The subsequent work added different variations to this model, several of them were reviewed and integrated in Stephan et al. (2004) [18] and Buxton (2004) [19]. Based on fundamental physiology, rather than empirical approaches, these enhanced models are able to unify existing literature and provide insight into how the underlying physiological mechanisms result in stable or/and transit BOLD response. However, the original model proposed by Buxton et al. and completed by Friston et al. is sufficient to account for the nonlinear behaviors observed in real-time series [5]. Too many state variables and parameters will not serve our purpose here better.

The dynamic intertwinement between multiple physiological variables, the cerebral blood flow (CBF) 𝑓, blood venous volume 𝑣, and veins deoxyghemoglobin content 𝑞, can be given as a set of nonlinear nondimensional differential equations [1, 20]: 𝑓̈=𝜖𝑢(𝑡)−𝑓̇𝜏_𝑠−𝑓−1𝜏_𝑓𝑣̇=1,𝑞̇=1𝜏₀𝑓−𝑣^1/𝛼,𝜏₀𝑓1−1−𝐸⁰^1/𝑓𝐸₀−𝑣^1/𝛼𝑞𝑣,(1) where 𝜖 is neuronal efficacy, reflecting the significance of neuronal activity evoked by experimental event, hence it varies with trial event; 𝜏_𝑠 and 𝜏_𝑓 represent time constant for signal decay and autoregulatory feedback from blood flow, respectively. The existence of feedback term can be inferred from the poststimulus undershoots in CBF [21]. The degree of nonlinearity of the BOLD signal is largely determined by the stiffness parameter 𝛼, which characterizes the balloon-like capacity of the venous compartment to expel blood at a greater rate when distended [22]. 𝐸₀ is resting net oxygen extraction fraction. All variables are expressed in normalized form, relative to resting values.

Noticing that the first equation has a second-order time derivative, so we can write this input-state-output system as a set of first-order ordinary differential equations by introducing another variable 𝑠=𝑓̇. By defining the state vector as 𝐱(𝑡)=[𝑓,𝑠,𝑞,𝑣]^𝑇, the system dynamic equation can be constructed from (1): 𝑥̇=𝑓(𝐱,𝜃,𝐮,𝐯)𝑣∼𝑁(0,𝐐),(2) where 𝜃={𝜖,𝜏_𝑠,𝜏_𝑓,𝜏₀,𝛼,𝐸₀,𝑉₀}∈ℝ^𝑙 is system parameter, the neuronal input 𝐮 represents system input, and 𝐯 is to account for the process noise.

The observed signal can be taken as a nonlinear function of volume 𝑣 and deoxyghemoglobin 𝑞 that comprises a volume-weighted sum of intra- and extravascular signal: 𝑦(𝑡)=𝑉₀𝑘¹(1−𝑞)+𝑘₂1−𝑞𝑣+𝑘³𝑘(1−𝑣)₁,=7𝐸₀,𝑘₂=2,𝑘₃=2𝐸₀−0.2,(3) appropriate for a 1.5 tesla magnet [1]. 𝑉₀ is the resting blood volume fraction, which generally varies across brain regions and subjects. All parameters are independent of each other. Their physiological definitions and probability distributions are given in Table 1 [23].

Table 1 

Hemodynamic model parameters and their probability distribution.

The actual observation 𝐲 is then composed of a deterministic part ℎ(𝑥,𝜃,𝑡) and a stochastic part 𝐰:𝐲=ℎ(𝐱,𝛽,𝐰)𝑤∼𝑁(0,𝐑),(4) where 𝐲 is the observation vector, 𝐰 is measurement noise, and 𝛽 consists of 𝑘₁,𝑘₂,and𝑘₃. Simultaneous estimation of 𝑉₀ and other parameters would be impossible, since their product is settled for each sampled measurement 𝑦_𝑘. The stiffness parameter 𝛼 is a nominal factor to BOLD contrast, it can be fixed to any value with its reasonable range in system identification [17].

Equation (2) describes a continuous-time hemodynamic process, and (4) models fMRI measurement as discrete sampling of the continuous system states, together they have formed a standard state-space representation for fMRI data assimilation. Given 𝑦_𝑘, the physiological states 𝐱 and the optimal parameters for a certain voxel can be estimated by use of UKF—a recursive minimum mean-square-error (MMSE) estimator.

2.2. Dual UKF and Joint UKF

The unscented Kalman filter has been applied extensively to the field of nonlinear estimation for both states and parameters. The basic framework of UKF involves estimation of the states of a discrete-time nonlinear dynamic system:𝐱𝐲_𝑘+1𝑘=𝐅𝐱_𝑘,𝐮_𝑘,𝐯_𝑘,=𝐇𝐱_𝑘,𝐧_𝑘,(5) where 𝐱_𝑘 represents the unknown system states, the system is driven by a known exogenous input 𝐮_𝑘 and process noise 𝐯_𝑘. The observation noise is given by 𝐧_𝑘.

The UKF generally involves recursive utilization of a deterministic “sampling’’ approach. The sampled points (sigma points) completely capture the true mean and covariance of the variables, and when propagated through the nonlinear system (𝐅 in this case), they are able to capture the posterior mean and covariance accurately to the 2nd order of Taylor series expansion [16].

Parameter estimation, or machine learning, on the other hand, involves determining a nonlinear mapping:𝐲_𝑘=𝐆𝐱_𝑘,𝐰,(6) where 𝐱_𝑘 is the input, 𝐲_𝑘 is the output, and the nonlinear map 𝐆(⋅) is parameterized by the vector 𝐰. Typically, a training set is provided by sample pairs consisting of known input and desired output, {𝐱_𝑘,𝐝_𝑘}. The goal of the learning can be expressed to some degree as solving for 𝐰 which minimizes the error of the machine: 𝐞_𝑘=𝐝_𝑘−𝐆(𝐱_𝑘,𝐰). In order to estimate the parameters by utilizing UKF, a new state-space representation can be written: 𝐝_𝑘𝐰_𝑘+1=𝐰_𝑘+𝐫_𝑘,=𝐆𝐱_𝑘,𝐰_𝑘+𝐞_𝑘,(7) where the parameters 𝐰_𝑘 correspond to a stationary process with identity state transition matrix, driven by process noise 𝐫_𝑘.

Given that BOLD signal is the only output and observation of the system in terms of Balloon model, the dual estimation problem, in which the system states and model parameters must be estimated simultaneously, can be given as follows:𝐱_𝑘+1𝐰=𝐅_𝑘+1𝐱_𝑘,𝐮_𝑘,𝐯_𝑘,𝐰_𝑘,𝐲_𝑘=𝐰_𝑘+𝐫_𝑘,=𝐆𝐱_𝑘,𝐰_𝑘+𝐞_𝑘.(8) Since standard UKF cannot be applied to this system immediately, dual UKF and joint UKF have been proposed as two alternatives. In the dual filtering method, two UKFs—one for state estimation, the other for parameter estimation—run in an alternate way. At each time step, the current estimate of parameters 𝐰_𝑘 is used in the state filter as given input, and likewise the current states estimate 𝐱̂_𝑘 is used in the parameter filter. On the contrary, a single UKF runs for both state and parameter estimation in the joint filtering. A higher dimensional joint state vector is defined: 𝐱̃_𝑘=[𝐱_𝑘𝑇𝐰_𝑘𝑇]^𝑇, and the state-space model is reformed as follows:𝐱̂𝐲_𝑘+1𝑘=𝐅𝐱̂_𝑘,𝐮_𝑘,𝐯_𝑘̂,=𝐆𝐱̂_𝑘,𝐧_𝑘.(9) The dual UKF and joint UKF approaches are illustrated in Figure 1.

Figure 1 

Schematic diagrams of joint filter (left) and dual filter (right).

In this section, the framework of UKF is briefly reviewed, a dual estimation problem with two approaches have been presented. In the next section, we will focus on examining the different performances of dual UKF and joint UKF, and all of our discussion will be in the context of Balloon model.

3. Results and Discussion

3.1. Biophysical Interpretation

One of the most prominent bifurcations between dual estimate and joint estimate is whether to incorporate interaction or correlation between states and parameters into filtering. As discussed earlier, the joint filter concatenates the state and parameter random variables into a single augmented state (𝑥̂̃_𝑘), so the cross-covariance between states and parameters is effectively modeled, that is, 𝐸𝑥̂̃−𝐸𝑥̂̃𝑥̂̃−𝐸𝑥̂̃^𝑇=⎡⎢⎢⎣𝐏^{𝑥𝑘𝑥𝑘}𝐏_{𝑥𝑘𝑤𝑘}𝐏_{𝑤𝑘𝑥𝑘}𝐏_{𝑤𝑘𝑤𝑘}⎤⎥⎥⎦.(10)

Dual filtering, on the other hand, decouples (in a statistical sense) the dual estimation problem by treating states and parameters separately, which means 𝐏_{𝑥𝑘𝑤𝑘}=𝐏_{𝑤𝑘𝑥𝑘}=0. For states and parameters involved in Balloon model, no dynamic interaction or biophysical correlation between them has been observed (they are uncorrelated variables), therefore, it is reasonable to expect dual filtering to exhibit more biophysical accuracy. Thus the fMRI experiments substantiated our assumption.

The real fMRI data was acquired from 8 health subjects. 136 acquisitions in total were made (RT = 2s), in block of 8, giving 16 16-second blocks. The condition for successive blocks alternated between rest and auditory stimulation, starting with rest. Auditory stimulation was emotionally neutral words presented at a rate of 60 per minute. We selected the largest activated voxels in superior temporal gyrus (GT) to implement data assimilation [24] (Figure 2). Bias correction was performed using the method in [25]. The two algorithms were initialized in identical way on experimental data and parameters.

Figure 2 

The greatest activated area of the group in the superior temporal gyrus (GT) for data assimilation.

Figure 3 shows the hemodynamic states given by joint UKF and dual UKF. The lower peak of blood flow inferred from joint UKF corresponds to the smaller neuronal efficacy (𝜖) in Figure 4.

(a) Joint estimate

(b) Dual estimate

Figure 3 

States estimated by dual and joint UKF. The dotted line corresponds to change of blood flow, the solid line shows venous volume and dashed line depicts dHb content.

(a) Joint estimate

(b) Dual estimate

Figure 4 

Parameters estimated by joint UKF and dual UKF. The mean values in Table 1 are used as initial values in our simulation. Apart from 𝐸0, which grows sharply at the beginning and does not change afterwards, all the parameters obtained from dual UKF can be seen as constants. Resting oxygen extraction 𝑉0 is the most important parameter in driving the model uncertainty [17], but simultaneous estimation of 𝑉0 and other parameters would be impossible, as stated earlier. 𝛼 is a nominal mechanism to BOLD signal. Therefore 𝑉0 and 𝛼 does not enter filtering.

Parameters estimated are shown in Figure 4. Signal decay, autoregulation et al. remain unchanged (almost) during dual filtering. While for joint UKF, the parameters do not converge to their final values until the 4th∼5th block (60∼80s after the first stimulation). This phenomenon is a strong indicator for the introduced interaction between states and parameters.

Real and estimated fMRI signals are plotted in Figure 5. The simulated BOLD signal given by dual UKF shows a slight overshoot, followed by gradual return to reduced plateau, and ending with a strong poststimulus undershoot. On the other hand, joint UKF fails to some degree to reconstruct a clean BOLD signal: the plateau is missing, neither does the evolving pattern of the signal show itself in a stable way in each block. The overestimated transit time (𝜏₀) leads to a reduction in amplitude of the BOLD peak; the less intense poststimulus undershoot can be explained by the underestimated signal decay (𝜏_𝑠). Given the fact that dual UKF and joint UKF have very similar performance for state estimation, which is made clear in Figure 3, it is safe to attribute the failure to the undesired fluctuations of parameters (especially 𝐸₀, which affects the estimated signal directly by (3)), or more precisely, the undesired interaction between parameters and states.

Figure 5 

Estimated and real fMRI signal.

3.2. Computational Interpretation

One of the most computationally expensive operations in UKF corresponds to calculating the new set of sigma points at each time update. This requires taking a matrix squareroot of the state covariance matrix, 𝐏∈ℝ^𝐿×𝐿, given by 𝐒𝐒^𝑇=𝐏 [16]. An efficient implementation using a Cholesky factorization requires in general 𝒪(𝐿³/6) computations [26]. Therefore enlarging the dimension of state-space vector will dramatically increase the computational complexity (also can be referred to as time complexity). In this subsection we will introduce two criteria for evaluating the property of dual UKF and joint UKF in time complexity.

Number of Floating-Point Operations (flops). In computing, floating point can be thought of as a computer realization of scientific notation, which is able to represent a wide range of values. Flops number required by a given algorithm or computer program is independent of the computing platform, although its precise value may differ under different counting rules. MATLAB (version before 6.0) has provided us with a useful function 𝑓𝑙𝑜𝑝𝑠 to specify the cumulative number of flops. For instance, if 𝐴 and 𝐵 are 𝑁-by-𝑁 matrixes, then the output of 𝑓𝑙𝑜𝑝𝑠(𝐴+𝐵) and 𝑓𝑙𝑜𝑝𝑠(𝐴𝐵) will approximately be 𝑁² and 2𝑁³.

Since Cholesky factorization is the only operation within UKF whose time complexity is proportional to 𝑁³, it is appropriate to consider that the flops number of Cholesky factorization is sufficient to determine the flops of the whole algorithm. For Balloon model, the dual estimation problem is about determining four state variables (𝐿_𝑥=4) and five parameters (all the parameters except 𝑉₀ and 𝛼, 𝐿_𝑤=5) at each time step. Table 2 shows the total flops number of Cholesky factorization involved in each predict-update cycle of UKF.

Table 2 

The fact that joint filtering requires half of the iterations that are required by dual filtering has been taken into account.

However, in practice by slightly restructuring the state vector, the process and observation models, we may introduce the noise with the same order of accuracy as the uncertainty in the state. First, the state vector is augmented to give a 𝐿^𝛼=𝐿_𝑥+𝐿_𝑣+𝐿_𝑛 dimensional vector;𝐱^𝛼𝑘=⎡⎢⎢⎢⎢⎣𝐱𝐯^𝑘_𝑘𝐧_𝑘⎤⎥⎥⎥⎥⎦.(11) Then the process model is rewritten as a function of 𝐱^𝛼𝑘, 𝐱^𝛼𝑘+1=𝐅(𝐱^𝛼𝑘,𝐮_𝑘); the unscented transform uses sigma points that are drawn from 𝐏^𝛼=⎡⎢⎢⎢⎢⎣𝐏^𝑥0𝐑0_𝑣000𝐑0_𝑛⎤⎥⎥⎥⎥⎦,(12) where 𝐑_𝑣 and 𝐑_𝑛 are the process and observation noise covariance. In this situation, similarly we can derive Table 3.

Table 3 

Flops number for augmented state vector.

For either case mentioned above, the flops number for joint filtering is at least 56% larger than that for dual filtering.

Comparing to flops count, overall execution time is a more tangible and practical criterion. We have tested our programs on several computers and collected their execution time data. Normally, a difference over 90% can be observed (16s and 30s for dual UKF and joint UKF, resp.). This result is even more impressive than that related to flops analysis, indicating that flops number is not the only factor influencing the operation time. Even if we take the rapid improvements in processing speed and memory into consideration, this variance is significant and should not be ignored.

4. Conclusions

In this paper we brought forward the dual Kalman filter as a reliable and efficient approach to estimating the states and parameters involved in balloon model. Comparing to the commonly used joint Kalman filter, its principle is in better conformity with the physiological reality, and by decoupling the dual estimation problem, it is much more calculational efficient to implement. The result of experiments showed good agreement with our conclusion.

Acknowledgments

This work is supported by the National Basic Research Program of China (no. 2010𝐶𝐵732500), the National Natural Science Foundation of China (no. 30800250), Doctoral Fund of Ministry of Education of China (no. 200803351022), Zhejiang Provincial Natural Science Foundation of China (no. 𝑌2080281), and Zhejiang Provincial Qianjiang Talent Plan (no. 2009𝑅10042).