Nonlinear estimation of BOLD signals with the aid of cerebral blood volume imaging

Functional magnetic resonance imaging (fMRI) offers a noninvasive technology to examine hemodynamic signals in the cerebrovascular system. The hemodynamic balloon model was introduced in 1998 to reveal the coupling dynamics between neural activity and blood oxygen level dependent (BOLD) responses by Buxton et al. [1]. The balloon model describes the causal mechanisms within a hemodynamic process in a certain region of interest (ROI) during brain activation. BOLD responses can be observed via the dynamic changes in cerebral blood volume (CBV) v, cerebral blood flow f, and vein deoxyhemoglobin (dHb) content q. This model is especially helpful in understand the potential consequences of interactions between physiological mechanisms. Since the inception of this model, there has been growing interest in using it to interpret observed fMRI data. The model can be used to infer biologically meaningful parameters that could be employed to investigate the changes in underlying physiological variables during brain activation [2–5], restrict the activation detection process with classic statistical techniques [6, 7], and deduce similar systems or different driving conditions [8–11].

The primary causes of unreliability in model estimation is that the BOLD fMRI technique is sensitive to changes in the signal from venous blood. The change in the signal intensity of a particular voxel is strongly dependent on what fraction of the voxel the vessel occupies. Moreove, changes in BOLD signal intensities during task activation are related not only to multiple physiological states but also regional vessel occupancy, including capillaries and large veins. Indeed, the evaluation of model structure also indicates that the blood volume fraction (BVF) greatly influences the uncertainty of model output [12]. However, this problem has been ignored in all previous studies. Most studies performed to data have avoided the ill-conditioning problem simply by employing a physiological plausible value of \(V_0=0.02\) instead of investigating the actual value in a particular ROI [2–5, 7, 13, 14] or throughout the brain [6, 15].

Given the importance of the true BVF, efforts are needed to incorporate actual vascular information of the voxel in the hemodynamic model estimation. Firstly, when a voxel includes only brain tissue, the assumption of \(V_0=0.02\) is reasonable [2, 16]. However, when a voxel is mostly or totally occupied by a vessel or vessels, the value might typically be above 0.6 [17]. Secondly, voxels associated with a larger amount of blood are always more likely to show significant BOLD activation due to the inherent nature of the fMRI technique. In this situation, employing an unrealistic \(V_0\) value might yield an unreliable model that does not reflect the physiological reality. This illustrates the importance of taking into account the actual BVF during the estimation procedure.

Several methods have been applied in attempts to obtain the true BVF. We recently showed that magnetic resonance angiography (MRA) might provide a method for roughly estimating the BVF value [18]. The results inferred that the \(V_0\) value in a voxel consists of two derivative components: (1) a constant tissue blood volume component \(V_s=0.02\), which is the small-vessel blood volume that includes capillaries and small postscapulaes, and (2) a variable large blood vessels component \(V_l\), which is the blood volume of large blood vessels. However, this method has not been used to obtain the actual \(V_0\) directly. Indeed, the regional CBV can be measured by another imaging modality, called the dynamic susceptibility contrast (DSC) material-enhanced gradient-echo (GE) MR technique [19]. The present study therefore augmented the true BVF acquired from CBV imaging in order to focuses on the influence of \(V_0\) on hemodynamic model estimation and the importance of using the true BVF in the analysis.

This paper is organized as follows. Firstly, we briefly review the hemodynamic Balloon model which constitutes the fundamental component of hemodynamic model estimation. Secondly, we explain the important influence of \(V_0\) with the adoption of a realistic value obtained from the CBV imaging technique. Lastly, the influence of \(V_0\) on model estimation within a single-region and multiple-regions according to the results of a classic bimanual finger tapping experiment is discussed in terms of the impacts of the actual \(V_0\) on parameter estimates and state-space reconstruction.

The hemodynamic balloon model describes the dynamic interrelationship between the blood flow f (neural activity to changes in flow), the regional blood volume information v (changes in flow to changes in blood volume and venous outflow), and the vein dHb content q (changes in flow, volume and oxygen extraction fraction to changes in dHb). The hemodynamic process can be described as the follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \ddot{f}=\epsilon u(t)-\frac{\dot{f}}{\tau _s}-\frac{f-1}{\tau _f} \\ \dot{v}=\frac{1}{\tau _0}(f-v^{1/\alpha }) \\ \dot{q}=\frac{1}{\tau _0}\left( f\frac{1-{(1-E_0)^{1/f}}}{E_0} -v^{1/{\alpha }}\frac{q}{v}\right) , \end{array}\right. } \end{aligned}$$

(1)

where \(\tau _s\) reflects signal decay, \(\tau _f\) is the feedback autoregulation time constant, \(\tau _0\) is the transit time, \(\alpha\) is a stiffness parameter, \(\epsilon\) is the neuronal efficacy, u(t) is the neuronal input, and \(E_0\) represents the resting oxygen extraction fraction. The variables f, v, and q are expressed in normalized form, relative to resting values. The balloon model can account for the hemodynamic responses in sparse, noisy fMRI measurements [12, 15]. However, since the above describing equations contain a second-order time derivative, we can introduce a new variable \(s=\dot{f}\) to express this hemodynamic system as a set of four first-order ordinary differential equations. Then the observed response BOLD signal could be expressed as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} y(t)=V_0(k_1(1-q)+k_2(1-\frac{q}{v})+k_3(1-v)), \\ k_1=7E_0, \quad k_2=2, \quad k_3=2E_0-0.2. \end{array}\right. } \end{aligned}$$

(2)

This equation is appropriate when using an fMRI machine 1.5-T magnet. The observed y is normalized relative to the value at rest, and \(V_0\) is the resting BVF [2]. Equations 1 and 2 consist of the architecture of hemodynamic input-and-output system. The model architecture is depicted in Fig. 1.

The BOLD response is associated with all of these parameters, but, we know that parameter \(V_0\) can not be identified along with other parameters simultaneously, instead only their product. Most previous efforts have imposed a physiologically plausible value of \(V_0= 0.02\) to handle the ill-conditioned nature of the problem [2–10]. Changes in the BOLD signal are strongly affected by \(V_0\), and so an unrealistic \(V_0\) may lead to unreliable model parameter estimation.

Two human subjects participated in this study. The experiment was approved by the Health Sciences Research Ethics Committee of Zhejiang University, and written informed consent was obtained from both subjects. Functional images were acquired on a 1.5-T scanner using a standard fMRI echo planar imaging protocol (resolution: \(64\times 64\) matrix; repetition time \(\mathrm{TR}=2\) s). In total, 110 acquisitions were made in a block-designed finger tapping experiment, giving 11 20-s blocks. The conditions for successive blocks alternated between rest and task performance, starting with rest. Furthermore, the CBV imaging sequence consisted of 30 T\(2^*\)-weighted images that were collected with a GE sequence (resolution: \(128\times 128\) matrix; 0.1 mmol/kg Gd-DTPA administered using a powered injector). In order to achieve a sufficient signal-to-noise ratio and complete coverage of the brain, TR was increased to 3.1 s, since a typical value is 1 s. The other sequence parameters remained unchanged.

All CBV images were down-sampled to make their spatial resolution identical to that of the fMRI image, and thereby allow voxel-by-voxel curve analysis. Concentration–time curves were created for each voxel [20–23]. The calculated \(V_0\) was then used in an existing data estimation procedure [24]. Figure 2 shows an example of an axial CBV image and the observed S and fitted concentration–time curves from one voxel. Data preprocessing and statistical analysis were performed using the SPM5 program (Wellcome Department of Cognitive Neurology, http://www.fil.ion.ucl.ac.uk/spm). The activation map was obtained by applying t-tests between all bimanual motor conditions and resting baselines with a cutoff for statistical significance of \(P < 0.001\).

This study focused on the important but long ignored issue of how the resting cerebral BVF (i.e., \(V_0\)) impacts hemodynamic models. Previous studies have used a physiologically plausible value of \(V_0= 0.02\) instead of exploring the actual \(V_0\) in the model estimation procedure. However, the intensity of any hemodynamic signal change is greatly affected by the regional BVF, since the active domains subject to model estimation often overlap with those areas characterized by a large BVF [30]. Under such circumstances, an inaccurate \(V_0\) may give rise to inaccurate estimates of the parameters and the reconstructed physiological state. This study used CBV imaging to augment the true \(V_0\) calculated in the hemodynamic model. In order to show the significance of applying the true \(V_0\), we have presented the different results obtained when using the real \(V_0\) and assumed \(V_0\) in terms of single-region model estimation and DCM. It was found that using the actual \(V_0\), yielded more realistic and physiologically meaningful model estimation results.

The results obtained in this study indicate that \(V_0\) has a rather complicated impact on estimated model parameters. Despite the BOLD responses being similar when using the assumed and real \(V_0\), there was a huge difference in the estimated parameters and the derived physiological state in the ROI. Because the balloon model describes the causal mechanism of a hemodynamic system, its order is higher than the externally observable system, which results in poorly identifiable model parameters due to the nature of nonlinear optimization and temporally sparse sampling. These model parameters have clear physiological meanings, and they should be justified and interpreted with caution [13, 31]. If the actual \(V_0\) is adopted, \(\epsilon\) can be more reliably observe via fMRI measurements. Therefore, \(V_0\) significantly influences the evaluations of brain connectivity. There have recently been extensive discussions on DCM and Granger causal modeling (GCM), with an emphasis on the connectivity among distributed brain systems [32–34]. In order to obtain a more robust understanding of brain causality, we used a biophysical model to eliminate signal bias in imaging procedure and variations of the hemodynamic response in diverse brain domains. However, an unrealistic \(V_0\) might degraded such efforts.

A potential limitation of the present study is to the extent that \(V_0\) as measured by CBV imaging is affected by the amount of blood associated with BOLD signals. We consider that both CBV imaging and the BOLD contrast have tiny difference in terms of the \(V_0\). The former contains the volume of blood across arteries, capillaries, and veins, whereas the latter is relevant to capillaries and veins [35]. Although the arterial fraction of CBV is much less than the venous BVF [36, 37], CBV imaging also partly removes the effect of overestimates about BVF. This is therefore a suitable method for approximating the value of \(V_0\). In addition, this study concentrated on explaining the influence of BVF on hemodynamic model estimation, and the results demonstrated the importance of taking advantage of actual BVF information in the estimation procedure. The argument about the origin of the two modalities were beyond the scope of this paper.

The present study presented the first empiric attempt to derive the actual \(V_0\) from data obtained using CBV imaging, with the aim of augmenting the existing estimation schemes. The results show that \(V_0\) significantly influences the estimation results within a single-region model estimation and DCM. Using the actual \(V_0\) can provide more reliable and accurate parameterizations and model predictions, and improve brain connectivity estimation based on fMRI data.