Exploiting Magnetic Resonance Angiography Imaging Improves Model Estimation of BOLD Signal

The change of BOLD signal relies heavily upon the resting blood volume fraction () associated with regional vasculature. However, existing hemodynamic data assimilation studies pretermit such concern. They simply assign the value in a physiologically plausible range to get over ill-conditioning of the assimilation problem and fail to explore actual . Such performance might lead to unreliable model estimation. In this work, we present the first exploration of the influence of on fMRI data assimilation, where actual within a given cortical area was calibrated by an MR angiography experiment and then was augmented into the assimilation scheme. We have investigated the impact of on single-region data assimilation and multi-region data assimilation (dynamic cause modeling, DCM) in a classical flashing checkerboard experiment. Results show that the employment of an assumed in fMRI data assimilation is only suitable for fMRI signal reconstruction and activation detection grounded on this signal, and not suitable for estimation of unobserved states and effective connectivity study. We thereby argue that introducing physically realistic in the assimilation process may provide more reliable estimation of physiological information, which contributes to a better understanding of the underlying hemodynamic processes. Such an effort is valuable and should be well appreciated.


Introduction
In 1998, Buxton and his colleagues introduced their celebrated hemodynamic model, Balloon model [1]. The comprehensive biophysical model of hemodynamic modulation describes the coupling dynamics from neural activity to observed blood oxygen level dependent (BOLD) signal [1,2]. It comprises the coupling mechanism of manifold physiological variables, blood flow (f ), blood volume (v), and deoxyhemoglobin content (q), during brain activation. This model then has been extended to include the effects of external inputs on blood flow inducing signal by Friston et al [3]. Since its inception, there is a growing interest in assimilating such a model with given sets of fMRI measurements in order to infer physiological parameters and associated states [4][5][6][7][8][9], constrain the activation detection process with classic statistics techniques [10,11], and extrapolate to similar systems and/or different driving conditions [12][13][14][15][16]. Although these works greatly enhance our understanding of the neural systems that mediate specific cognitive processes, they are still kind of problematic in offering reliable inference on the hemodynamic system behaviors.
The query on reliability of estimation primarily comes from the assumption about resting blood volume fraction (V 0 ) in the assimilation procedure. It has long been noted that BOLD contrast is highly weighted by venous blood content. The change of signal intensity in given region thereby depends heavily on local vessel geometry including capillaries and large veins. The evaluation of model structure also indicates that V 0 is a leading influence mechanism in driving the model output uncertainty [17]. However, This parameter can not be identified along with other model parameters simultaneously due to the ill-conditioning of the inverse problem. All studies so far have engaged a physiological plausible value V 0~0 :02 in region of interest (ROI) [3][4][5][6][7][8][9]18] or throughout the whole brain [10,11] to dispel the ill-conditioning problem, instead of investigating actual V 0 . When a voxel includes only brain tissue, the assumption V 0~0 :02 is reasonable [3,19]. When a voxel is mostly or totally occupied by a vessel or vessels, however, the value might typically be above 0:6 [20]. On the other hand, these voxels that contain large blood content are always more likely to show significant BOLD activation due to the nature of fMRI technique. In this situation, the employment of unrealistic V 0 value in data assimilation might produce unreliable model estimation, far straying from physiological reality. The effort of incorporating actual vascular information of voxels into the fMRI data assimilation therefore should be well appreciated.
In this study, we presented the first attempt to exploit actual resting blood volume fraction in assimilation procedure. The actual V 0 is derived from the segmentation of the vein in the MR angiography (MRA), then augmented into the existing assimilation schemes. As physical realistic V 0 is adopted in assimilation process, more reasonable inference about hemodynamic behavior can thus be expected. We will illustrate the efficacy of the combinative approach on single-region data assimilation and multi-region data assimilation.
This paper is organized as follows. We first simply review the hemodynamic Balloon model and its formulation that forms the basis of data assimilation, then describe the derivation of V 0 from MRA images. The impacts of actual V 0 on states forecast and parameter estimation are presented in terms of data assimilation and dynamic causal models subsequently.

Hemodynamic Balloon Model
The original hemodynamic Balloon model consists of three subsystem linkings: (1) neural activity to changes in flow; (2) changes in flow to changes in blood volume and venous outflow; (3) changes in flow, volume and oxygen extraction fraction to changes in deoxyhemoglobin (dHb). It describes the dynamic intertwinement between the blood flow f , the blood venous volume v and the veins dHb content q, which can be given as the following [1,3]: where E is neuronal efficacy; u(t) is the neuronal input; t s reflects signal decay; t f is the feedback autoregulation time constant; t 0 is the transit time; a is the stiffness parameter; and E 0 denotes the resting oxygen extraction fraction. All variables are expressed in normalized form, relative to resting values. The input-state-output system is represented by nonlinear equations of a series of physiological states. Equation (1) has a second-order time derivative, and we can write this system as a set of four first-order ordinary differential equations (ODEs) by introducing a new variable s~_ f f . Although the Balloon model is enhanced somehow afterwards [21][22][23], the model structure analysis shows that the original model is sufficient to account for the hemodynamic response in sparse, noisy fMRI measurement [11,17].
Furthermore, the BOLD observation can be expressed as: appropriate for a 1.5 Tesla magnet [1]. V 0 is the resting blood volume fraction, which may vary across brain regions and across subjects. The model architecture is summarised in Figure 1.
For any given combination of parameters q~fE,t s ,t f , t 0 ,a,E 0 ,V 0 g and neuronal inputs u, equations (1) and (2) produce a predicted BOLD response. They form the basis for fMRI data assimilation from the measured dataset. Note that parameter V 0 can not be identified along with other parameters simultaneously, but only their product is admitted. Up to now, all existing efforts have circumvented the ill-conditioning nature by imposing a physiological plausible value V 0~0 :02 [4][5][6][7][8][10][11][12][13][14][15][16]. However, the usage of this parameter value is expedient so that the assimilation problem can be solved. Since the change in the BOLD signal depends heavily on V 0 , unrealistic V 0 may lead to unreliable model parameter estimation. In addition, the stiffness parameter a shows a marginal influence to the system output variance, and it can be fixed within its physiological reasonable range (a~0:33 here) without significant loss of information in data assimilation processing [8,10,11,17].

Derivation of V 0 from MR Angiography Image
In hemodynamic model, V 0 is defined as the venous volume of blood present in a voxel. It represents the ratio occupied vessels with sizes ranging from capillaries to large veins that all contribute to fMRI measurements in the area [1,2]. Typical resting value in brain tissue which only contains capillaries is around 2 per cent. When a vessel or vessels are present in a voxel, local blood volume will dramatically rise. The value in large vessel region is typically above 0:6. The presence of large vessel is expected to make V 0 inhomogeneous. Fortunately, large blood vessel are accessible by MR angiography (MRA) imaging.
Consider that V 0 in a voxel consists of two different derivative components, constant tissue blood volume component V s and varied large blood vessels component V l : V s~0 :02 is small-vessel blood volume (including capillaries and small postcapillaries). V l is blood volume of large blood vessels (veins and venules). It is associated with draining veins, and spatially varies across different brain areas in general. In this study we made use of high-resolution time-of-flight magnetic resonance angiography (TOF-MRA) scanning to accurately locate the blood vessels in the brain. The principle of TOF-MRA imaging is based on the enhancement of the signal of dynamic blood flow and the suppression of the signal of static tissues. The resolution of the TOF-MRA image was 0:9375|0:9375|1:5mm 3 (intensity range [0,1425]), which was much better than that of the fMRI image (3:75|3:75|6mm 3 , in this study). All images were collected in the same field of view (FOV). In TOF, veins usually bear higher signal intensities than the surrounding tissues, thus making the segmentation of major veins feasible and reliable. It is practicable to downsample the fine vasculature information to coarse fMRI scale in order to obtain the estimation of regional V l at given fMRI voxels. We therefore attempted to combine the MR angiography image and the fMRI image for uncertainty reduction in data assimilation. SPM2 program (Wellcome Department of Cognitive Neurology, http://www.fil.ion.ucl.ac.uk/spm) was used for our data preprocessing, voxel by voxel. Each fMRI volume was realigned to the first volume, and created a mean of the realigned data. The mean functional image then was upsampled to the resolution of the MRA. The MRA image was coregistered (Estimation and Reslice in SPM software) to the resultant upsampled mean image with linear interpolation. For high SNR TOF images, we performed the segmentation by simple thresholding. In the experiment, the segmentation threshold was set to 300 by simple visual guidance. After the vascular segmentation, we can obtain the large blood vessel composition (i.e. V l ) of each voxel in EPI images. Moreover, an isolated voxel with intensity higher than 300 was considered as noise and was therefore excluded from the calculation of V 0 . Since the MRA image has a much higher spatial resolution than fMRI image (64 : 1 in this study), the large vessel fraction V l of each voxel was expressed as: Combined with the small-vessel fraction V s , the total blood volume fraction of each voxel in fMRI image was expressed: The first term represents the volume of blood from large vessels in a voxel, the second term is the blood volume fraction in the remaining brain tissue. In this sense, MRA can be thought of as an indirect, physical measurement of V 0 , and can be treated as 'true' V 0 value. V 0~0 :02 is a special case of the formulation when large vessel does not exist in this voxel.

Experiment and Data Preprocessing
The participant provided written informed consent before beginning the experiment which was approved by the Health Sciences Research Ethics Committee of Zhejiang University. Functional images were acquired on a 1:5-Tesla scanner (Marconi EDGE ECLIPSE) using a standard fMRI gradient echo echo-planar imaging (EPI) protocol (TE, 40 ms; TR, 2000 ms; flip angle, 90 0 ; NEX, 1; FOV, 24 cm; resolution, 64|64 matrix). Sixteen contiguous 6-mm-thick slices, 0:5-mm-intervals were acquired to provide a coverage of the entire brain. Foam padding was used to limit head motion within the coil.
Block design experiment was performed in this study. The subject was presented with classical flashing checkerboard pattern when scans were acquired. Activation maps (Pv0:01) were generated with the SPM software package (Wellcome Department of Cognitive Neurology, http://www.fil.ion.ucl.ac.uk/spm), which used a General Linear Model approach to detect regions with significant response during the task.
The Time-Of-Flight MR Angiography (TOF) image was segmented to extract the major veins in the brain. In TOF, veins usually bear higher signal intensities than the surrounding tissues, which makes the segmentation of major veins easy. For the high signal-noise-ratio (SNR) TOF images, we performed the segmentation by simple thresholding. The selection of threshold could be accomplished manually. In our work, we segmented veins by thresholding because of the high quality of the TOF image. Figure 2 presents an example of CBV imaging segmented from one subject. After thresholding, a de-noise step such as an opening operation, could be executed in order to eliminate some isolated noises. Once the segmented vasculature was obtained, we need to further transfer the information of vein position to the fMRI data. This requires that we registered/aligned the TOF to the fMRI image. Since the subject of the two images was the same patient in a short period, multi-modal rigid registration was enough to perform the task. We chose the classic mutual information as the similarity metric to do the 3D registration in the physical domain.
After finishing the segmentation and registration, we obtained corresponding brain blood volume in the voxel. The actual V 0 was then augmented into existing assimilation schemes [8,11,24]. In this study, two ROIs were selected from the greatest activation locus of primary visual cortex (V1) and V4 (Figure 3). We defined the clusters based on edges, but not included corners, so that the voxel had 4 neighbors (in the same slice). Figure 3 clearly shows

Results
The Impact of V 0 on Single-region Data Assimilation In this study, as a demonstration, we chose V1 as region of interest. We estimated the state functions and model parameters ( Figure 3). Since there are not large veins in V4 area, this approach makes no difference in this area. For the sake of simplicity, we assumed a constant neural parameter E 1~E2~Á Á Á~E n through-out all trials, where n denotes trial number. The estimation scheme is formally identical to that in [24][25][26]. Figure 4 shows reconstructed BOLD response (left) and underlying physiological states (right) in the greatest activation locus of primary visual cortex (Figure 3) with actual V 0 , given as solid line. As a comparison, we also evaluated the estimated BOLD response and physiological states with assumed V 0~0 :02 value, which was widely employed in previous studies [3][4][5][6][7][8][10][11][12][13][14][15][16]18], given as dash line in Figure 4. We found that two different V 0 values produced very similar BOLD estimates (Figure 4, left), only tiny discrepancy in post stimulus undershoot stage could be found. Nevertheless, a significant distinction was observed in reconstructed physiological states between two values (Figure 4, right). Though the experimental stimulus induced a puny change in the blood flow f , the blood venous volume v and the veins dHb content q, the approach used assumed V 0 deduced a substantial  change during task due to magnifying effect of large blood content. This implied that the presence of large veins in an activated area contributed excess signal in this area. The change of BOLD signal in this area mainly derived from the large-vessel signal, not from the multiple physiological states, namely, not from the experimental related neuronal activity. Since statistical inference essentially is grounded on the amplitude of BOLD response, this area may surely be considered active in statistic analysis of the BOLD signal change, though it is absent at the response efficacy elicited by neuronal activity. This explains why the employment of assumed V 0 in detection process still could generate very similar activation map with those obtained from classic linear model [10,11]. The same difference also can be found in estimated model parameters (Table 1). Specially, neuronal efficacy E is 0:0094 with actual V 0 derived from MRA image, while E is 0:1534 with assumed V 0 value. Assumed, underestimated V 0 substantially overestimates the neuronal efficacy parameter, E. Since parameter E reflects the efficacy with which neuronal activity causes an increase in signal, we argue that the estimated efficiency parameter E in each voxel might be a good index to sign actual activation level.
The Impact of V 0 on Dynamic Causal Models (DCM) Dynamic causal modeling (DCM) has been introduced as a generic method to explore effective connectivity from the hemodynamic observations [12,27]. Apart from Balloon model, this model additionally embeds a neurobiological modelization of the dynamic interactions among brain areas into the hemodynamic models in these areas, and it can be regarded as an extension of hemodynamic model from single region to covering multiple regions. Single-region data assimilation supposes that extrinsic experimental input consistently accesses all brain regions, whereas DCM designs that inputs produce responses in two different ways: extrinsic influence from sensory input and intrinsic influence from interaction regions. As uncertain V 0 makes the greatest impact on estimates of neuronal efficacy parameter E in hemodynamic model, it is interesting to investigate the effect of V 0 on DCM.
In this study, as an example, two regions were selected using maxima of activation map to construct the hierarchical system. The system architecture was shown in Figure 5. The two maxima were located in visual area V1 and V4. Region-specific time series comprised all neighbor voxels of each maxima location (a total of 5 voxels). The location is shown in Figure 3. The system describes a simple hierarchy of forward connections where two primary motor regions influence each other, and can be expressed as the following where x 1 (x 2 ) is neuronal dynamics in V1 (V4); u 1 and u 2 represent external inputs into the system; a 11 (a 22 ) represents the inner connectivity within the region in the absence of input; a 12 (a 21 ) encodes the fixed inter-regions connectivity in the absence of input; c 11 (c 22 ) embodies the extrinsic influences of input on neuronal activity. Equations (6) and (7) then were appended into the states vector [8,11]. The measurement vector was expanded to include two observations in the two areas as well. Two inputs corresponded to a 0{1 quarewave function for the occurrence of experimental stimulus ( Figure 5). The outputs of the system are two time series from two regions. The estimation scheme employed for DCM is formally identical to that in previous studies [8,11].
The results of this analysis are presented in Figure 5. The connections are shown as directed black arrows with the coupling parameters calculated with actual V 0 alongside. The values in brackets are parameters estimated with assumed V 0 . V 0~0 :4947 in visual area V1, V 0~0 :02 in V4 and assumed V 0~0 :02 in two areas. As expected, the significant difference in connectivity parameters with actual V 0 and assumed V 0 can be found ( Figure 5). The fixed connectivity from V1(V4) to V4(V1) is {0:28 while considering the contribution of vessels, whereas the value is {0:50 while the effect was discounted.
From the above analysis, the employment of an assumed V 0 in the hemodynamic data assimilation seems to be only suitable for fMRI signal reconstruction and activation detection grounded on this estimated signal, not for effective connectivity study that by means of estimated neuronal activity (e.g. Eu(t)) makes inference about the coupling among brain areas and how that coupling is influenced by changes in experimental context. Due to the regulation of resting blood volume fraction (V 0 ), in fMRI imaging, large BOLD signal changes are often associated with large draining veins, while tissue areas have low BOLD signal changes. These results suggest that the impact of V 0 on fMRI data assimilation should be considered. Actual V 0 should be investigated or these areas that are dominated by large veins should be excluded in the region-specific analysis.

Discussion
This work is principally concerned with an important but long ignored issue in previous efforts on the hemodynamic model -the influence of resting cerebral blood volume fraction V 0 . Previous studies postulated a physiologically plausible value V 0~0 :02 in assimilation procedure to handle ill-posedness of the problem, as opposed to exploring true BVF. This practice may lead to inaccurate results. In this study, instead of arbitrary assignment, we propose a combinative approach that supplements realistic V 0 derived from MR angiography (MRA) image into an existing hemodynamic assimilation scheme to achieve more reliable model estimation. We find that V 0 has a complicated influence on assimilation results. Though arbitrarily assigned V 0 can produce similar BOLD response with realistic V 0 , there is significant difference in reconstructed physiological states and estimated model parameters, indicating that the application of these parameter estimated by assumed V 0 should be justified and interpreted with caution [7,28]. Moreover, as uncertain V 0 value leads to larger deviation in estimated efficacy parameters E than that in other parameters, we also have investigated the influence of V 0 on dynamic causality modeling which estimates connectivity in different brain regions by means of E estimates. Not surprisingly, Table 1. Estimated model parameters with true value (V t ) and typical assumed value (V a ) in the greatest activated locus of primary visual cortex (V1). V 0 has also considerable impact on the evaluation of brain connectivity. We thereby argue that introducing more realistic V 0 in DCM can provide more reliable estimation of interregional coupling, and assist to acquire a better understanding of brain connectivity that is of considerable interest in neuroimage community recently, such as Human Connectome Project (HCP) in NIH, Brain CONNECT Project in Europe, and National Basic Research Program of China (973) under Grant 2011CB707800.
A possible criticism on this work is to what extent MRA image is able to provide accurate actual cerebral blood volume fraction reflecting BOLD response. Indeed MRA imaging is not a direct, physical measurement of V 0 . However, as noted, the MRA image has a much higher spatial resolution than a fMRI image (64 : 1 in this study). In this sense, MRA can be thought of as an indirect, physical measurement of large veins component V l in given areas. Combining with tissue blood volume component V s , a more realistic V 0 value can be obtained. An imperfect measurement is better than arbitrary assignment without any measurement. As more physical realistic V 0 is incorporated into the assimilation procedure, more reliable information of the underlying physiological dynamics is reconstructed, and physiologically more meaningful results may be expected. Another limitation is that few experiments were performed in this study. The main cause is the incorrect but pervasive belief that MRI scans are harmful to health in China. The recruitment of subjects was very difficult. This was even reported recently by Nature [29]. However, our study is mainly concerned with the influence of V 0 on assimilation, and to discuss the importance of introducing actual V 0 information into the assimilation process. The results clearly illustrate our intention. We therefore believe that more experiments is not necessary. Despite these limitations, we argue that such a effort is valuable and should be well appreciated, in particular, while nearly 400 studies on hemodynamic data assimilation have been reported every year [30]. Currently, we are trying to experimentally verify that augmenting more realistic V 0 derived from MRA imaging into assimilation process will provide more accurate states forecast and parameter estimation.