Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2012 (2012), Article ID 961967, 7 pages
Research Article

Reliable and Efficient Approach of BOLD Signal with Dual Kalman Filtering

State Key Laboratory of Modern Optical Instrumentation, Zhejiang University, Hangzhou 310027, China

Received 25 May 2012; Accepted 11 July 2012

Academic Editor: Huafeng Liu

Copyright © 2012 Cong Liu and Zhenghui Hu. 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.


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. [1113] 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 𝒪(𝐿3) [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𝑣=𝜏0𝑓𝑣1/𝛼,1̇𝑞=𝜏0𝑓11𝐸01/𝑓𝐸0𝑣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]. 𝐸0 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 𝜃={𝜖,𝜏𝑠,𝜏𝑓,𝜏0,𝛼,𝐸0,𝑉0}𝑙 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: 𝑦(𝑡)=𝑉0𝑘1(1𝑞)+𝑘2𝑞1𝑣+𝑘3,𝑘(1𝑣)1=7𝐸0,𝑘2=2,𝑘3=2𝐸00.2,(3) appropriate for a 1.5 tesla magnet [1]. 𝑉0 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 𝑘1,𝑘2,and𝑘3. Simultaneous estimation of 𝑉0 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.

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.
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 4th5th block (6080s 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 (𝜏0) 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 𝐸0, 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 𝒪(𝐿3/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 𝑁2 and 2𝑁3.

Since Cholesky factorization is the only operation within UKF whose time complexity is proportional to 𝑁3, 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 𝑉0 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 𝐏𝛼=𝐏𝑥000𝐑𝑣000𝐑𝑛,(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.


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).


  1. R. B. Buxton, E. C. Wong, and L. R. Frank, “Dynamics of blood flow and oxygenation changes during brain activation: the balloon model,” Magnetic Resonance in Medicine, vol. 39, no. 6, pp. 855–864, 1998. View at Publisher · View at Google Scholar · View at Scopus
  2. R. D. Hoge, J. Atkinson, B. Gill, G. R. Crelier, S. Marrett, and G. Bruce Pike, “Investigation of bold signal dependence on cerebral blood flow and oxygen consumption: the deoxyhemoglobin dilution model,” Magnetic Resonance in Medicine, vol. 42, pp. 849–863, 1999. View at Google Scholar
  3. Y. Zheng, J. Martindale, D. Johnston, M. Jones, J. Berwick, and J. Mayhew, “A model of the hemodynamic response and oxygen delivery to brain,” NeuroImage, vol. 16, no. 3, pp. 617–637, 2002. View at Publisher · View at Google Scholar · View at Scopus
  4. T. Obata, T. T. Liu, K. L. Miller et al., “Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: application of the balloon model to the interpretation of BOLD transients,” NeuroImage, vol. 21, no. 1, pp. 144–153, 2004. View at Publisher · View at Google Scholar · View at Scopus
  5. K. J. Friston, A. Mechelli, R. Turner, and C. J. Price, “Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics,” NeuroImage, vol. 12, no. 4, pp. 466–477, 2000. View at Publisher · View at Google Scholar · View at Scopus
  6. K. J. Friston, “Bayesian estimation of dynamical systems: an application to fMRI,” NeuroImage, vol. 16, no. 2, pp. 513–530, 2002. View at Publisher · View at Google Scholar · View at Scopus
  7. T. Deneux and O. Faugeras, “Using nonlinear models in fMRI data analysis: model selection and activation detection,” NeuroImage, vol. 32, no. 4, pp. 1669–1689, 2006. View at Publisher · View at Google Scholar · View at Scopus
  8. J. J. Riera, J. Watanabe, I. Kazuki et al., “A state-space model of the hemodynamic approach: nonlinear filtering of BOLD signals,” NeuroImage, vol. 21, no. 2, pp. 547–567, 2004. View at Publisher · View at Google Scholar · View at Scopus
  9. L. A. Johnston, E. Duff, and G. F. Egan, “Particle filtering for nonlinear BOLD signal analysis,” Lecture Notes in Computer Science, vol. 4191, pp. 292–299, 2006. View at Google Scholar · View at Scopus
  10. L. A. Johnston, E. Duff, I. Mareels, and G. F. Egan, “Nonlinear estimation of the BOLD signal,” NeuroImage, vol. 40, no. 2, pp. 504–514, 2008. View at Publisher · View at Google Scholar · View at Scopus
  11. Z. H. Hu and P. C. Shi, “Nonlinear analysis of bold signal: biophysical modeling, physiological states, and functional activation,” in Proceedings of the 10th Internation Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '07), pp. 734–741, Brisbane, Australia, 2007.
  12. P. Shi, Z. Hu, X. Zhao, and H. Liu, “Nonlinear analysis of the BOLD signal,” Eurasip Journal on Advances in Signal Processing, vol. 2009, Article ID 215409, 2009. View at Publisher · View at Google Scholar · View at Scopus
  13. Z. Hu, C. Liu, P. Shi, and H. Liu, “Exploiting magnetic resonance angiography imaging improves model estimation of bold signal,” PloS One, vol. 7, no. 2, Article ID e31612, 2012. View at Google Scholar
  14. K. J. Friston, N. Trujillo-Barreto, and J. Daunizeau, “DEM: a variational treatment of dynamic systems,” NeuroImage, vol. 41, no. 3, pp. 849–885, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. J. Daunizeau, K. J. Friston, and S. J. Kiebel, “Variational bayesian identification and prediction of stochastic nonlinear dynamic causal models,” Physica D, vol. 238, no. 21, pp. 2089–2118, 2009. View at Publisher · View at Google Scholar · View at Scopus
  16. R. van der Merwe and E. A. Wan, “The square-root unscented Kalman filter for state and parameter-estimation,” in Proceedings of the IEEE Interntional Conference on Acoustics, Speech, and Signal Processing (ICASSP '01), vol. 6, pp. 3461–3464, May 2001. View at Scopus
  17. Z. H. Hu and P. C. Shi, “Sensitivity analysis for biomedical models,” IEEE Transactions on Medical Imaging, vol. 29, no. 11, pp. 1870–1881, 2010. View at Publisher · View at Google Scholar · View at Scopus
  18. K. E. Stephan, L. M. Harrison, W. D. Penny, and K. J. Friston, “Biophysical models of fMRI responses,” Current Opinion in Neurobiology, vol. 14, no. 5, pp. 629–635, 2004. View at Publisher · View at Google Scholar · View at Scopus
  19. K. Uludaǧ, R. B. Buxton, D. J. Dubowitz, and T. T. Liu, “Modeling the hemodynamic response to brain activation,” NeuroImage, vol. 23, supplement 1, pp. S220–S233, 2004. View at Publisher · View at Google Scholar · View at Scopus
  20. R. B. Buxton and L. R. Frank, “A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation,” Journal of Cerebral Blood Flow & Metabolism, vol. 17, no. 1, pp. 64–72, 1997. View at Google Scholar · View at Scopus
  21. K. Irikura, K. I. Maynard, and M. A. Moskowitz, “Importance of nitric oxide synthase inhibition to the attenuated vascular responses induced by topical L-nitroarginine during vibrissal stimulation,” Journal of Cerebral Blood Flow & Metabolism, vol. 14, no. 1, pp. 45–48, 1994. View at Google Scholar · View at Scopus
  22. D. E. Glaser, K. J. Friston, A. Mechelli, R. Turner, and C. J. Price, “Haemodynamic modelling,” Human Brain Function, Elsevier Press, London, UK, 2003. View at Google Scholar
  23. R. L. Bellaire, E. W. Kamen, and S. M. Zabin, “A new nonlinear interated filter with application to target tracking,” in Proceedings of the AeroSense: 8th International Symposium on Aerospace/Defense Sensing, Simulation, and Controls, vol. 2561, pp. 240–251, 1995.
  24. Z. Hu, P. Ni, C. Liu, X. Zhao, H. Liu, and P. Shi, “Quantitative evaluation of activation state in functional brain imaging,” Brain Topography. In press. View at Publisher · View at Google Scholar
  25. Z. Hu, H. Liu, and P. Shi, “Concurrent bias correction in hemodynamic data assimilation,” Medical Image Analysis. In press. View at Publisher · View at Google Scholar
  26. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 2 edition, 1992.