Research Article  Open Access
A Sensitivity Analysis of fMRI Balloon Model
Abstract
Functional magnetic resonance imaging (fMRI) allows the mapping of the brain activation through measurements of the Blood Oxygenation Level Dependent (BOLD) contrast. The characterization of the pathway from the input stimulus to the output BOLD signal requires the selection of an adequate hemodynamic model and the satisfaction of some specific conditions while conducting the experiment and calibrating the model. This paper, focuses on the identifiability of the Balloon hemodynamic model. By identifiability, we mean the ability to estimate accurately the model parameters given the input and the output measurement. Previous studies of the Balloon model have somehow added knowledge either by choosing prior distributions for the parameters, freezing some of them, or looking for the solution as a projection on a natural basis of some vector space. In these studies, the identification was generally assessed using eventrelated paradigms. This paper justifies the reasons behind the need of adding knowledge, choosing certain paradigms, and completing the few existing identifiability studies through a global sensitivity analysis of the Balloon model in the case of blocked design experiment.
1. Introduction
Functional Magnetic Resonance Imaging (fMRI) technique has been the area of great interest since it measures the activity of the brain, which gives much insight on brain functions. The contrast agent in fMRI is based on the rate of oxygen consumption in the brain through the Blood Oxygenation Level Dependent (BOLD) signal. Whether the aim is to characterize a region of the brain or to compare different regions and/or subjects, it is necessary to link the BOLD indirect measure to the stimulus.
Several studies have considered modeling the hemodynamic response of the brain. These models fall in two categories: general linear modelingbased detection methods, which assume a known shape of the hemodynamic response function (HRF) and estimation ones which track the most appropriate HRF form in regions eliciting evoked activity in response to stimulus presentation. The first class of models uses statistical tools (discrepancy) to fit some basis functions (Poisson function [1], Gaussian function [2], Gamma function [3], inverse Logit function [4]) to the BOLD for the aim of activation detection. The second is more concerned with the physiological aspects that underlie the BOLD transients. It is natural at this point to think that understanding the brain functions requires a much more elaborated framework to map the neural activity to changes in the measured BOLD signal than a convolution by some predefined hemodynamic response functions. Furthermore, it has been shown [5–7] that to be more consistent, this mapping has to encompass nonlinear effects to describe more accurately the pathway from neural activation to the BOLD.
Among the nonlinear hemodynamic models, linking changes in the physiological variables, namely, the blood flow, the volume, and the oxygenation level to the neural activity in a localized area of the brain, is the Balloon model. This model has been first proposed by Buxton et al. [6]. Several variants of the Balloon model have been developed such as Friston et al.’s [8] variant considered in this paper. The main feature of the Balloon model is that it describes the BOLD signal as a nonlinear combination of the blood volume and deoxyhemoglobin concentration, which are expressed as nonlinear functions of the neural activity and the blood flow.
Models in fMRI, in particular the Balloon model, have been used in many studies to estimate underlying parameters. Indeed, these models become criticized because (1) our comprehension of the mechanisms underlying the hemodynamic response is evolving and (2) they are often used without previous evaluation of their identifiability. Therefore, this paper is relevant to this second point. Indeed, estimating the state, the input, and parameters for the Balloon model using BOLD measurements is illposed. The illposedness, also seen as an identifiability issue when the input signal is given, consists in having a small change in the BOLD signal while observing large variations of parameters and thus poor accuracy of the parameter estimates. Different strategies have been used to add prior knowledge that makes the problem wellposed. Riera et al. [9] used radial basis functions to express the input while Havlicek et al. [10] supposed Gaussian priors for the parameters. Small prior variances are sought to overcome the illposedness of some parameters. In the same field of Kalman filtering techniques, Hu et al. [11] used a square root unscented Kalman (SRUKF) filter to estimate the states and some parameters with Gaussian priors.
Besides, even if it were possible to develop theoretical tools to get an analytical expression for the parameters from the nonlinear equations describing the dynamic behavior of the Balloon model, the input signal has to be “rich enough” so that the parameters can be estimated accurately. The typical example for that is when one parameter quantifies the output behavior with respect to the derivative of the input, and if the input is constant, then the parameter cannot be estimated efficiently. Thus, the input sequence has to be well designed to capture the effect of each of the unknown parameters on the output signal.
In practice, usually there is unfortunately a tradeoff between the good contrasttonoise ratio (due to slow or static behavior) and the estimation efficiency (often assessed for fast time scales). In general, the challenge is attributed to experimental limitations such as possible damage of the studied system (due to high switchings in frequencies, e.g.) and poor signals quality (low signal to noise ratio SNR) when dealing with high frequencies. This question has been raised for the fMRI experimental paradigms by Liu et al. [12] and a theoretical model for the quantification of the socalled detection power/estimation efficiency duality. The approach was developed within the hemodynamic general linear modeling framework.
In this study, we will focus on the identifiability study of the Balloon model parameters which represents an attempt to answer the question: Given the output measurement (BOLD) and the input (Neural activity), can (Balloon model parameters) be determined uniquely?
The paper is organized as follows. In the second section, we introduce the nonlinear hemodynamic Balloon model and present the commonly used fMRI paradigms. The third section reviews the existing sensitivity approaches of the Balloon model and completes them by a global sensitivity study when the stimulus follows block design experiment. Finally some conclusions are presented in the fourth section.
2. Problem Statement
2.1. The Balloon Model
The dependence of the BOLD signal upon the variations of the normalized blood flow , volume , and deoxyhemoglobin content has been proposed by Buxton et al. [6] and completed by Mandeville et al. [13]. The model presents the venous compartment’s structure as a Balloon where a stimulus resulting in a local neuronal activity leads to an increase in the blood flow as illustrated in Figure 1.
This increase exceeds the cerebral metabolic rate of oxygen thus reducing the concentration of deoxyhemoglobin which in turn results in an increase in the magnetic resonance signal, also called the BOLD effect [14]. Although the proposed model provides a mapping between the BOLD and the physiological variables, there still remains the need to explicit the relationship between the neural activity and the blood flow dynamics. The modeling of this part of the hemodynamic mechanism has been introduced by Friston et al. [8]. Indeed, under the assumption of a linear relation between the synaptic activity and the regional cerebral blood flow, a new variable called the flowinducing signal has been introduced. The proposed coupling between the neuronal activity and the cerebral blood flow has a linear second order behavior of a relaxed oscillator.
To sum up, the Balloon model can be formulated as a single input/single output nonlinear statespace representation:where the vector is a state vector which consists of a set of the nonmeasurable variables defined above and is called the “flowinducing signal.” The physical hidden states (, , and ) are normalized with respect to the rest phase and thus equal to at rest. is the input of the model which represents the stimulus function, generally associated to the neural activity, and is the output given by the BOLD signal. and are nonlinear evolution and observation functions which are given by represents a set of parameters and consists of the neural efficiency , the flow decay , the flow time constant , the venous transit time , Grub’s stiffness parameter , the oxygen extraction at rest , and the blood volume fraction at rest . These parameters may vary from a brain region to another or from a subject to another. The parameters , , and depend on the scanner.
The main feature of the Balloon model is that it is nonlinear both in state and parameters. It presents an unusual nonlinearity with respect to the state vector component, for and dynamic equations.
2.2. Experimental Design in fMRI
Experiments in fMRI are sought to test biological hypotheses. Nevertheless, there are several experimental constraints related to both physiological effects to which fMRI is sensitive and temporal properties of the BOLD signal itself [15]. The resulting noises include, for example, the cardiac and respiratory cycles as well as motion artifacts. For a more thorough analysis about the nature of fMRI signals’ components, readers can refer, for example, to Howseman et al. [16]. These limitations restrict the frequency of stimuli occurrence that provides an unambiguously detectable activation signal. According to the mentioned study, the duration of a cycle ranges from few minutes to few seconds. The cycles consist of successive boxcar functions where the rest periods are intended to provide a constant baseline.
2.2.1. Block Paradigm
In this category [17, 18], the stimulus consists of alternated periods of activation and rest. Each task is of roughly equal duration (the typical block length is 15 s on and 15 s off). In this case, the activation signal and the noise correlation are low and a large BOLD signal change [19] is introduced. The rest periods provide lowfrequency drift secondary to physiological noise which can be considered as baseline.
2.2.2. EventRelated Paradigm
This experimental design, requiring a higher temporal resolution, emerged in the mid 90s taking advantage of faster image acquisition techniques. The event durations are much shorter (about one or two seconds) [18] compared to block design paradigm. The main advantages of eventrelated design include the reduction of the effect of head movements and the ability to capture transients in the BOLD signal [19]. The interstimulus length can be as short as the event duration or may be varied to avoid the subject predict the events and more importantly to improve the detectability of the effects of interest.
2.2.3. Mixed Paradigm
This paradigm results from the combination of block and eventrelated designs. It consists of close events separated by long rest periods. Mixing transient (trialrelated) and sustained (taskrelated) components in the obtained BOLD signals allows the identification of regions of the brain with moderelated neural activity and others with taskrelated patterns.
3. Identifiability Study of the Balloon Model
Deneux and Faugeras [20] studied the identifiability of the Balloon model parameters at a local level by addressing, in a neighborhood of some nominal values for the parameters , the invertibility of the mapping: where is the measured BOLD signal, is the neural activity, considered as a known input, and represents the set of unknown parameters.
The rank of the Jacobian of with respect to determines the number of “unknowns” that can be estimated accurately. These unknowns can match the parameters if the directions defined by the Jacobian are orthogonal; otherwise a small change in one parameter can be partially (if the Jacobian is full rank) or completely offset by a change of a combination of the remaining parameters. Thus to quantify the sensitivity of the output to a small change in a parameter , an orthogonal “projection” of the Jacobian with respect to on that of is considered. In other words, ifthenwhere is the pseudo inverse of , which is taken as a measure of the sensitivity of to a small change in .
Although this sensitivity study is at a local level and it is predominately oriented to eventrelated like input , the obtained results for two sets of neural activity signals (a single impulse and two successive impulses) are interesting. They show that and estimations are more accurate while seems to be the most difficult to tackle.
More recently, Hu and Shi [21] led to a global sensitivity analysis on some biomedical systems including the Balloon model. This method, developed in a Bayesian framework, consists in assigning a prior distribution to each of the model parameters and then quantifying its contribution to the model output through the study of the posterior distribution. Given the measurement and an unknown variable , then the first order contribution of to the output is quantified as follows: where is the expectation operator and the variance operator. Note that the sum of all is always less or equal to 1; the equality is assessed in the case of additive models. Higher order contributions represent joint variance terms with the remaining variables. For example, second order contributions, also called second order sensitivity indexes are terms of the form:and so forth. According to this criterion the higher the index is the more likely the variable to drive the model output variance is, and thus the better the candidate for model calibration is. Similarly, variables with small indexes can, for example, be frozen to their nominal values. In practice, this criterion requires evaluation of probability densities for nonlinear (nonadditive) quantities which makes this technique computationally costly to implement. Two approximations based on Monte Carlo method and the Unscented Transformation were used.
Sensitivity analysis has been applied to the Balloon model in the case of eventrelated and block designs. The main outcomes of this study show that (the Balloon model parameter, not to be confused with the above notations for variances) has the highest importance in driving the BOLD uncertainty, whereas seems to have a minor effect making it a noninfluential parameter for the BOLD signal. Regarding the dynamics of the BOLD signal, seems to affect more the initial dip and the rising to peak phase while contributes more at the undershoot level.
In order to come up with a conclusion on the identifiability of the Balloon model and to consolidate the previous results for the special case of block design experiments, the next part of this work will be dedicated to a sensitivity study.
3.1. Variation of the BOLD Signal When Varying One Parameter
This is an intuitive approach which consists in tracking the changes in the output signal, that is, the BOLD when only one of the model’s parameters varies. It allows seeing to which parameters the BOLD is the more sensitive. For simplicity, we will name the Balloon model parameters in (2) with a change to the inverse for and . This parametrization together with the nominal values taken for the reference model (as in Friston et al. 2000 [8]) are summarized in Table 1.

As for the remaining observation equation parameters, it is agreed that they are scanner dependent [6]. The commonly used values, corresponding to , are taken as in the literature:
Remark. and are dependent on , so they represent additional unknowns to the estimation problem.
The input signal consists of a square function with an amplitude of one and a period of 1 minute (30 seconds OFF and 30 seconds ON). The sampling period is equal to s. The input and the BOLD signal are shown in Figure 2. The experiment corresponds to short transients compared to the overall stimulus duration and the BOLD signal presents static regime intervals which is a typical feature of block designs.
Each subplot of Figure 3 (except the first one which is the nominal BOLD) corresponds to the relative error in percent of the BOLD signal when one of the parameters takes values in a grid around the nominal value (shown in Table 1).
The Minimum. All the plots, except that of , present a unique extreme in a neighborhood of the nominal value of the parameter corresponding to the global minimum. Subsequently, the estimation of can be more complex, requiring a good minimum search method and an adequate initial guess.
Symmetry of the Criterion. It is also worth noticing that the criterion is not symmetric for , , , and . Indeed, the greater the value of the parameter is, the smaller its contribution to the criterion is. This makes it obvious that, to estimate these parameters efficiently, initial guesses should be taken less than the expected value for each of the parameters.
The dissymmetry of the cost function might be explained as follows: these parameters monitor the transients of the BOLD signal in the sense that small values imply slow dynamics for the Balloon model. A slow system is unable to track changes of the input which results in a big disparity between the nominal and actual BOLD signals. Similarly, a faster dynamic system reaches the steady state more quickly than the nominal system. The disparity only applies to the transients, which are very short compared to steady state regimes.
Sensitivity of the BOLD. The shape of the criterion suggests that the BOLD is less sensitive to changes in the abovementioned parameters when their values exceed the nominal values. As an illustration, a variation of more than of is needed to increase the relative error of the BOLD signal by . While the latter value of the error is a very satisfying rate when fitting the BOLD signal, it is highly inaccurate to estimate a parameter with uncertainty.
Apart from the contribution of each parameter to the BOLD, it is also interesting to check if there is a correlation between these single contributions. In this case, the Balloon model, subject to blocked paradigm, is defined by the different correlation directions rather than the parameters.
To “partially” answer this question, the contributions of parameters, two by two, to the BOLD signal are explored in the next subsection.
3.2. Sensitivity of the BOLD to Simultaneous Changes in Two Parameters
The approach adopted in this part of the study consists in creating a grid for each of the parameters. Each plot corresponds to the relative error of the BOLD, when a pair of parameters varies on the grid and the remaining parameters are maintained to their nominal values. The grids are chosen such that the relative error on the BOLD does not exceed , assuming that any efficient approach for the estimation of the model parameters will not provide a worse result for the BOLD fitting.
The pairs having less correlated contributions are summarized in Figures 4 and 5. It is interesting to notice that the lowest correlation is obtained for the neural efficiency () and the venous transit time () as shown in Figure 5(b). More generally, the venous transit time seems to define a direction that is independent from contributions of most of the Balloon parameters (more precisely, it only slightly correlates with ).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
If the global minimum is well marked for small variations around the nominal values (corresponding to a relative error less than ), correlation directions are more and more clearly defined as the relative error increases. These correlations are almost exclusively located at regions corresponding to parameter values which are greater than the nominal value, which is typically the case for () and (). This further confirms that the BOLD signal is less sensitive to large increase in the parameters defining the dynamics of the Balloon model.
The second series of plots, shown in Figures 6 to 8, illustrates cases where the global minimum is not “sufficiently marked” even in a small neighborhood of the nominal values of parameters. Clearly, the estimation of is problematic. Indeed, whatever parameter contribution, taken in a reasonable neighborhood of the nominal value, can be offset by .
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusion
The sensitivity study of the Balloon model highlights the importance that should be given to the design of the experimental paradigm. This work completes the previous studies, showing the difficulty to estimate accurately the Balloon model parameters using the wellknown blocked design. The lack of accuracy is due in part to the tradeoff between detectability and estimation when it comes to fMRI experiments. Indeed, increasing the randomness of the input signal, as usually required to cover as much frequencies as possible, makes it difficult to have a contrast signal which is high enough to detect the active areas of the brain.
The need to develop new paradigms that could mix up randomness and longer events seems obvious especially in the absence of a wellestablished theoretical framework that can handle the exact identifiability study of the Balloon model. The calibration could also be assessed through several experiments, eventrelated like inputs, for example, allow the estimation of dynamic parameters (while blocked designs allow the detection of evoked activity more reliably). The estimation of , however, remains an issue as it can easily introduce a significant bias on the overall estimates.
More generally speaking, this study goes beyond the Balloon model itself to further confirm the difficulty of characterizing the pathway from the stimulus to the BOLD effect. Several open questions, tightly connected to the inability of the BOLD signal to calibrate the hemodynamic model accurately, can be raised. Possible solutions include revising the hemodynamic modeling, or combining fMRI with other imaging techniques to have extra measurements, for example, the cerebral blood flow.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST). The authors would like to thank anonymous reviewers for valuables comments on the paper.
References
 K. J. Friston, P. Jezzard, and R. Turner, “Analysis of functional MRI timeseries,” Human Brain Mapping, vol. 1, no. 2, pp. 153–171, 1993. View at: Google Scholar
 J. C. Rajapakse, F. Kruggel, J. M. Maisog, and D. Y. von Cramon, “Modeling hemodynamic response for analysis of functional MRI timeseries,” Human Brain Mapping, vol. 6, no. 4, pp. 283–300, 1998. View at: Publisher Site  Google Scholar
 P. Ciuciu, J.B. Poline, G. Marrelec, J. Idier, C. Pallier, and H. Benali, “Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment,” IEEE Transactions on Medical Imaging, vol. 22, no. 10, pp. 1235–1251, 2003. View at: Publisher Site  Google Scholar
 M. Lindquist and T. Wager, “Modeling the hemodynamic response function using inverse logit functions,” in Proceedings of the Human Brain Mapping Annual Meeting, June 2005. View at: Google Scholar
 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 and Metabolism, vol. 17, no. 1, pp. 64–72, 1997. View at: Google Scholar
 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 Site  Google Scholar
 K. J. Friston, O. Josephs, G. Rees, and R. Turner, “Nonlinear eventrelated responses in fMRI,” Magnetic Resonance in Medicine, vol. 39, no. 1, pp. 41–52, 1998. View at: Publisher Site  Google Scholar
 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 Site  Google Scholar
 J. J. Riera, J. Watanabe, I. Kazuki et al., “A statespace model of the hemodynamic approach: nonlinear filtering of bold signals,” NeuroImage, vol. 21, no. 2, pp. 547–567, 2004. View at: Publisher Site  Google Scholar
 M. Havlicek, K. J. Friston, J. Jan, M. Brazdil, and V. D. Calhoun, “Dynamic modeling of neuronal responses in fMRI using cubature Kalman filtering,” NeuroImage, vol. 56, no. 4, pp. 2109–2128, 2011. View at: Publisher Site  Google Scholar
 Z. Hu, X. Zhao, H. Liu, and P. Shi, “Nonlinear analysis of the BOLD signal,” Eurasip Journal on Advances in Signal Processing, vol. 2009, Article ID 215409, 2009. View at: Publisher Site  Google Scholar
 T. T. Liu, L. R. Frank, E. C. Wong, and R. B. Buxton, “Detection power, estimation efficiency, and predictability in eventrelated fMRI,” NeuroImage, vol. 13, no. 4, pp. 759–773, 2001. View at: Publisher Site  Google Scholar
 J. B. Mandeville, J. J. A. Marota, C. Ayata et al., “Evidence of a cerebrovascular postarteriole windkessel with delayed compliance,” Journal of Cerebral Blood Flow and Metabolism, vol. 19, no. 6, pp. 679–689, 1999. View at: Google Scholar
 R. B. Buxton, K. Uludag, D. J. Dubowitz, and T. T. Liu, “Modeling the hemodynamic response to brain activation,” NeuroImage, vol. 23, no. 1, pp. S220–S233, 2004. View at: Publisher Site  Google Scholar
 R. B. Buxton, Introduction to Functional Magnetic Resonance Imaging Principles and Techniques, Cambridge University Press, Cambridge, UK, 2nd edition, 2009.
 A. M. Howseman, O. Josephs, G. Rees, and K. J. Friston, “Special issues in functional magnetic resonance imaging,” SPM Course, Short Course Notes, 1997. View at: Google Scholar
 A. M. Blamire, S. Ogawa, K. Ugurbil et al., “Dynamic mapping of the human visual cortex by highspeed magnetic resonance imaging,” Proceedings of the National Academy of Sciences of the United States of America, vol. 89, no. 22, pp. 11069–11073, 1992. View at: Publisher Site  Google Scholar
 R. L. Buckner, P. A. Bandettini, K. M. O'Craven et al., “Detection of cortical activation during averaged single trials of a cognitive task using functional magnetic resonance imaging,” Proceedings of the National Academy of Sciences of the United States of America, vol. 93, no. 25, pp. 14878–14883, 1996. View at: Publisher Site  Google Scholar
 E. Amaro Jr. and G. J. Barker, “Study design in fMRI: basic principles,” Brain and Cognition, vol. 60, no. 3, pp. 220–232, 2006. View at: Publisher Site  Google Scholar
 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 Site  Google Scholar
 Z. Hu and P. Shi, “Sensitivity analysis for biomedical models,” IEEE Transactions on Medical Imaging, vol. 29, no. 11, pp. 1870–1881, 2010. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Chadia Zayane and Taous Meriem LalegKirati. 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.