Volumetric BOLD fMRI simulation: from neurovascular coupling to multivoxel imaging
 Zikuan Chen^{1}Email author and
 Vince Calhoun^{1, 2}
DOI: 10.1186/14712342128
© Chen and Calhoun; licensee BioMed Central Ltd. 2012
Received: 28 December 2011
Accepted: 23 April 2012
Published: 23 April 2012
Abstract
Background
The blood oxygenationlevel dependent (BOLD) functional magnetic resonance imaging (fMRI) modality has been numerically simulated by calculating single voxel signals. However, the observation on single voxel signals cannot provide information regarding the spatial distribution of the signals. Specifically, a single BOLD voxel signal simulation cannot answer the fundamental question: is the magnetic resonance (MR) image a replica of its underling magnetic susceptibility source? In this paper, we address this problem by proposing a multivoxel volumetric BOLD fMRI simulation model and a susceptibility expression formula for linear neurovascular coupling process, that allow us to examine the BOLD fMRI procedure from neurovascular coupling to MR image formation.
Methods
Since MRI technology only senses the magnetism property, we represent a linear neurovascularcoupled BOLD state by a magnetic susceptibility expression formula, which accounts for the parameters of cortical vasculature, intravascular blood oxygenation level, and local neuroactivity. Upon the susceptibility expression of a BOLD state, we carry out volumetric BOLD fMRI simulation by calculating the fieldmap (established by susceptibility magnetization) and the complex multivoxel MR image (by intravoxel dephasing). Given the predefined susceptibility source and the calculated complex MR image, we compare the MR magnitude (phase, respectively) image with the predefined susceptibility source (the calculated fieldmap) by spatial correlation.
Results
The spatial correlation between the MR magnitude image and the magnetic susceptibility source is about 0.90 for the settings of T_{E} = 30 ms, B_{0} = 3 T, voxel size = 100 micron, vessel radius = 3 micron, and blood volume fraction = 2%. Using these parameters value, the spatial correlation between the MR phase image and the susceptibilityinduced fieldmap is close to 1.00.
Conclusion
Our simulation results show that the MR magnitude image is not an exact replica of the magnetic susceptibility source (spatial correlation ≈ 0.90), and that the MR phase image conforms closely with the susceptibilityinduced fieldmap (spatial correlation ≈ 1.00).
Keywords
Bold fMRI Neurovascular coupling Neuroactive blob (NAB) Blood magnetism Intravoxel dephasing Voxelization Magnetic fieldmap Spatial correlationBackground
Blood oxygenationlevel dependent (BOLD) functional magnetic resonance imaging (fMRI) has been widely accepted for brain functional mapping and neuroimaging [1–6]. The imaging principle of BOLD fMRI is that: a neuroactivity incurs cerebral vascular blood magnetism perturbation that can be detected by magnetic resonance imaging (MRI). In neurophysiology, the BOLD fMRI can be described by a neurovascular coupling model as follows [7–11]: a neuronal activity incurs a vascular response in terms of changes in cerebral blood flow (CBF), cerebral blood volume (CBV), and cerebral metabolic rate of oxygen (CMRO_{2}). The neuroactivityinduced biomagnetic susceptibility perturbation can be detected by T2*weighted MRI (T2*MRI) [12]. An overall BOLD fMRI model can be decomposed into a cascade of two modules [13]. One is a neurophysiology module that addresses the vascular response to a neuroactivity in context of neurovascular coupling [7–11]; The neurovascular coupling process produces an intravascular blood biomagnetism perturbation that is detectable by T2*MRI [12]. Another is a MRI technology module dedicated to imaging the susceptibilityexpressed neurovascular coupling states. Despite its wide acceptance for neuroimaging and brain mapping, the BOLD fMRI mechanism is not fully understood [5, 14, 15]. Specifically, the complicated neurovascularcoupled BOLD process is not fully understood, and the imaging performance of T2*MRI detection on a BOLD state has never been quantitatively examined in a multivoxel simulation.
In past decades, there have been many published reports on numerical simulations of BOLD mechanism [1, 6, 16–18]. To our knowledge, previous BOLD simulations were carried out for investigating single voxel signals and did not address multiple voxels. This may be partially due to the computational burden and the implementation difficulty in manipulating a large 3D image matrix (e.g. 2048 × 2048 × 2048 in our simulation). In practice, a BOLD fMRI experiment produces a 4D dataset that consists of a time series of 3D image matrices (typically in size of 64 × 64 × 32 voxels) captured by T2*MRI at a discrete time point (snapshot). A multivoxel image offers many spatial features, such as "contrast and texture", "geometry: edge/shape/pattern", "topography", "topology", and so on. Furthermore, given an input source and its corresponding output image, we can assess its overall imaging performance in terms of point spread function, spatial invariance/variance, and linearity/nonlinearity of a digital imaging system.
Since T2*MRI is designed to sense an inhomogeneous fieldmap that is established via magnetization of an inhomogeneous susceptibility distribution, the underlying source of fMRI is the susceptibilityexpressed distribution of a neurovascular coupling state. Therefore, for numerical BOLD fMRI simulation, we need numerically characterize the neurovascular coupling process in terms of biomagnetism susceptibility perturbation for the purpose of T2*MRI detectibility. Given a susceptibility map (representing a snapshot of dynamic BOLD susceptibility perturbation), we can carry out T2*MRI simulation in a way similar to BOLD voxel signal simulations [19–21]. Intuitively, a multivoxel BOLD fMRI simulation may be implemented by spatially arranging individual voxel signal values. In the implementation, nonetheless, we need account for the electromagnetic interaction among the voxels (due to the nonlocal effect of vascular blood magnetization [22]), especially for the magnetic influence from vasculatures in surrounding voxels that may influence the BOLD voxel signal by 10% [23].
In this paper, we propose a volumetric BOLD fMRI model that deals with a cortical field of view (FOV), as a whole matrix (without dividing into submatrices), in which the magnetic influence among the intraFOV vasculatures are accounted for during the fieldmap calculation using a Fourier technique [16, 23–25]. After calculation of the fieldmap over the FOV, we proceed to calculate the output image of T2*MRI based on intravoxel dephasing. Upon the forward BOLD fMRI simulation, we obtain a complexvalued image (consisting of a pair of magnitude and phase parts) at each snapshot time. By comparing the MR magnitude image with the susceptibility source in a measure of spatial correlation, thereby evaluating how well the susceptibility source can be reproduced in the MR magnitude image. Meanwhile, we may compare the MR phase image with the fieldmap, thereby justifying the practice of fieldmap measurement by MR phase imaging. It is noted that the fieldmap is spatially different from the susceptibility source by a 3D convolution [13, 22, 26], therefore, we do not expect a morphological match between the susceptibility map and either the fieldmap or the phasemap.
The neurovascular coupling process is a complicated neurophysiological process. Experiments have shown that the BOLD activity in response to a neuroactivity may be described as a linear response model [7, 8, 27–32]. For the implementation of numerical simulation, we propose a linear neurovascular coupling formula to account for the biomagnetic susceptibility contributions from various neurovascular coupling aspects, such as neuroactivity (NAB: neuroactivity blob), vasculature (CBV: cerebral blood volume), blood physiology (CBF: cerebral blood flow), and blood metabolism (CMRO_{2}:cerebral metabolism rate of oxygen).
A BOLD fMRI study produces a 4D MR dataset, of which each 3D volume at a time point is interpreted as a snapshot of a dynamic BOLD state. The 4D BOLD fMRI data acquisition can be numerically simulated by a train of 3D T2*MRI snapshots provided the dynamic susceptibility perturbation is numerically specified for each snapshot time. In this report, we will only focus on the snapshot capture of a susceptibilityexpressed BOLD physiological state at a specific time point. Our approach is general and we can implement a dynamic 4D BOLD fMRI simulation by repeating the snapshot imaging at a series of time points at more computation cost.
Methods
Neurovascular coupling formulation
where (σ_{ x }, σ_{ y }, σ_{ z } )delineates an ellipsoidal profile (σ_{ x } = σ_{ y } = σ_{ z } ) for sphere, σ_{x} = σ_{y} << σ_{z} for long ellipsoid or cylinder, σ_{x} ≠ σ_{y} ≠ σ_{z} for general ellipsoid), (x_{0}, y_{0}, z_{0}) denotes the location of the NAB in the FOV, and c = max(NAB) represents the maximum activity at the NAB center (the activity strength is scaled by 0 < c < 1). We should mention that a nonspherical blob like the one in Figure 1 can be formulated by a combination of several regular Gaussianshaped NAB primitives that are different in (x_{0}, y_{0}, z_{0}) and (σ_{x}, σ_{y}, σ_{z}). The graded local neuroactivity strength over the FOV is reflected in the spatially distributed multivalues in the range of [0, c] (maximum activity strength at the NAB center, moderate strength at the NAB boundary, and zero strength outside the NAB).
where dom(V) denotes the space domain of the selected cortical vasculature, and bfrac(t) the dynamic blood volume fraction, and V the volume of FOV. It is noted that the time parameter t is explicitly used to indicate the variation in vasculature geometry. We should point out that we can use bfrac(t) to numerically characterize the blood physiological parameters: CBV and CBF. Specifically, a vasodilation associated with CBV manifests as a slight increase in bfrac(t), that is bfrac(t+Δt) > bfrac(t) for Δt > 0; and an accelerated blood flow can be equivalent to an increase in bfrac(t) as well. Since bfrac plays a control parameter during random vasculature generation (in Eq. (2)), its effect on intravascular blood physiology is exerted through the vasculature configuration geometry.
It is also known that only the red blood cells in blood stream convey oxygen that contribute to intravascular blood susceptibility perturbation. The total volume of red blood cells in normal blood is about 40%, as described in terms of hematocrit (Hct≈0.4). Blood physiology also shows that a red blood cell can carry up to 4 oxygen molecules (via attachment to 4 heme groups in a hemoglobin). Due to the oxygen detachment during microcirculation, the deoxygenated blood reveals more paramagnetism than the oxygenated blood, that is χ_{deoxy} > χ_{oxy}. Let Y(t) represent the dynamic blood oxygenation level, then (1Y(t)) represents the dynamic deoxygenation level, that is a parameter to reflect the cerebral metabolism of oxygen (CMRO_{2}). It is noted that Y(t) = [0,1], with Y = 1 for the fully oxygenated blood in artery and Y = 0 for the fully deoxygenated blood in vein.
Parameters and settings for numerical simulations
Parameters  Settings  Remarks 

B_{0}  3 Tesla  Main static magnetic field 
Hct  0.4 (dimensionless ratio)  Blood hematocrit 
Y  0.6 (dimensionless ratio)  Oxygenation level = [0,1] 
χ_{do}  0.27 × 4π ppm (SI metrics)  = χ_{dexoy}  χ_{oxy} 
NAB  NAB(x, y, z) Gaussian ellipsoid  Neuroactive blob (Eq(1)) 
V(x, y, z)  Cortical vasculature (Eq(2))  
FOV: D_{0} × D_{0} × D_{0}  2 × 2 × 2 mm^{3}  Cortical field of view 
FOV support matrix  2048 × 2048 × 2048  Digital geometry (1 μm grid) 
voxel: d_{0}× d_{0}× d_{0}  32 × 32 × 32, 64 × 64 × 64, 128 × 128 × 128  Three voxel sizes (1 μm grid) 
C[x_{n}, y_{n}, z_{n};T_{E}]  64 × 64 × 64, 32 × 32 × 32, 16 × 16 × 16  Complex multivoxel image 
A[x_{n}, y_{n}, z_{n};T_{E}]  Assuming nonnegative  MR magnitude image(Eq(9)) 
P[x_{n}, y_{n}, z_{n};T_{E}]  Assuming positive/zero/negative  MR phase image (Eq(9)) 
corrA  In a range of [0, 1]  Magnitude mapping (Eq(10) 
corrP  In a range of [0, 1]  Phase mapping (Eq(11)) 
T_{E}  In a range of [0, 60] ms  Gradient echo time 
 1.
It is a linear neurovascular coupling formula that accounts for different neurovascular parameters by a spatiotemporal modulation. For simplicity, we do not consider the hemodynamic time lag, spatial displacement, or spatial response spread in this work, though our model will support it.
 2.
The BOLD susceptibility perturbation is due to the temporal modulation by the blood deoxygenation level (1Y(t)), which is an embodiment of CMRO_{2}. The interplay among CBF, CBV and vasculature are numerically characterized by a single parameter bfrac(t) that reflects the blood dynamic through vasculature configuration (V(x, y, z, t)).
 3.
Only intravascular blood deoxyhemoglobin contributes to the BOLD susceptibility perturbation; no contribution from extravascular tissue (due to V(x, y, z) = 0 for extravascular region), nor from oxyhemoglobin (due to 1Y = 0 for Y = 1), nor from neuronal inactive regions (due to NAB = 0).
 4.
The volumetric computational model can be considered as a generalization of the single voxel neurovascular coupling model (Δχ = χ _{do}⋅Hct⋅(1Y)) that has been accepted for single voxel BOLD signal simulation [1, 3, 17, 23].
Forward BOLD fMRI simulation
From magnetism perturbation to fieldmap establishment
where FT and IFT stand for Fourier transform and inverse Fourier transform respectively, (k_{x}, k_{y}, k_{z}) coordinates in the Fourier domain, conv the convolution, and h(x, y, z) the kernel of magnetic point dipole field.
Multivoxel partition of FOV
where * denotes the convolution, and V [x_{n}, y_{n}, z_{n}] connotes the voxelization of V(x, y, z). It is shown that the voxelization operation (denoted by squared brackets "[ ]") suppresses the intravoxel details and produces a digital image representation of the vasculature configuration in the FOV. By applying voxelization to the 3D distributions of NAB(x, y, z), Δχ(x, y, z), ΔB(x, y, z), we obtain the corresponding 3D matrices of NAB[x_{n}, y_{n}, z_{n}], Δχ[x_{n}, y_{n}, z_{n}], and ΔB[x_{n}, y_{n}, z_{n}], respectively. For example, in our simulation (see later), the fieldmap is originally represented as a large matrix in size of 2048 × 2048 × 2048 (with a fine grid resolution as used for the digital FOV representation), which can be reduced a smaller matrix in size of 32 × 32 × 32 by voxelization with voxel size of 64 × 64 × 64.
Voxel signal calculation by intravoxel dephasing integration
where the echo time T_{E} is explicitly retained to remind of the T_{E} dependence of voxel signal (as will be demonstrated in our simulation later).
From voxel signal values to a multivoxel image
Where C[x_{n}, y_{n}, z_{n}; T_{E} = 0]denotes the nondecay initial magnitude image and ∠C[x_{n}, y_{n}, z_{n}; T_{E} = 0] the initial phase image. Henceforth, we will denote the MR magnitude and phase images as A[x_{n}, y_{n}, z_{n}; T_{E}] and P[x_{n}, y_{n}, z_{n}; T_{E}], respectively.
Backward mappings
In the forward BOLD fMRI calculation, we obtain a pair of MR magnitude image and phase image by Eq. (9). Considering the NAB[x_{n}, y_{n}, z_{n}] (the voxelized version of the NAB(x, y, z) in Eq. (6)) as the neuronal origin, the conventional neuroimaging effort consists in establishing a backward mapping from the MR magnitude images to the neuronal origin, as designated by "A ~ NAB" in Figure 2. Under a linear neuruovascular coupling model, we simplify the mapping between the vascular response (Δχ[x, y, z]) and the neuronal origin (NAB[x, y, z]) by a linear mapping in Eq. (4). On one hand, Δχ represents the vascular response to the neuronal origin under a neurovascular coupling process; on the other hand, it plays a vascular origin for the T2*MRI detection. In this paper, we will focus on the backward mapping from the MR magnitude image to its vascular origin as designated by "A~Δχ" in Figure 2.
Where cov(x, y) denotes the covariance between vector x and y, std(x) the standard deviation, and ":" a nDto1D operator (as used by Matlab language) that reorders a highdimensional array entries into one dimensional vector. The correlation coefficient defined in Eq. (10) gives rise to corrA ∈ [0, 1] with corrA = 1 for a perfect match (linear mapping) and corrA ≠1 for mismatch (nonlinear mapping). Since there is little decay for a short relaxation time (A[x_{n}, y_{n}, z_{n}; T_{E}]→0 for T_{E} → 0), a meaningful corrA(T_{E}) should be evaluated at a relative long T_{E} (T_{E} > 0). However, a long T_{E} will also introduce a diffusion smearing effect, which is not addressed herein. It is noted that corrA defined in Eq. (10) is a numerical measure of the backward mapping "A~Δχ" in Figure 2, which connotes the magnitudevssusceptibility correlation.
Likely, corrP is a numerical measure of the phasevsfieldmap correlation, a backward mapping designated by "P~ΔB" in Figure 2.
In summary, our computational BOLD fMRI model can be used for backward mappings: magnitudevssusceptibility correlation (corrA) and phasevsfieldmap correlation (corrP) by Eqs. (10) and (11), respectively. The goal of conventional neuroimaging and brain mapping is to render a backward mapping from MR magnitude image to neuronal origin, which consist of two steps: from the MR magnitude image to its vascular origin (as designated by "A~Δχ" and numerically measured by corrA) and then from the vascular response to its neuronal origin. In this paper, we are concerned with the mapping of "A~Δχ", with which we show the effect of MRI technology on the imaging performance of BOLD fMRI.
Results
The main parameters and their settings for the simulations we performed are listed in table 1.
The overall numerical simulation scheme is described by a flowchart in Figure 2. It starts with defining a NAB in a cortical FOV (Eq. (1)) and filling the FOV with randomly generated vascular networks (Eq. (2)). For digital geometry depiction of small vessels, the FOV is finely gridded and represented in a large support matrix in size of 2048 × 2048 × 2048 (grid resolution = 1 micron). Then the magnetic susceptibility expression of a neurovascularcoupled BOLD state is analytically described by a linear spatiotemporal modulation formula in Eq. (4). At a snapshot time, we predefine a 3D susceptibility distribution based on the following parameter settings: a Gaussianshaped NAB, a vasculatureladen FOV (under control of bfrac(t) = 2%) and Y(t) = 0.6. Let B_{0} = 3 T, we calculate the fieldmap from a 3D susceptibility perturbation in size of 2048 × 2048 × 2048 (assuming the same support matrix as FOV) by Eq. (5). By spatially partitioning the FOV into voxels, we calculate the voxel signals for a T_{E} by using the intravoxel dephasing integration in Eq. (7). At last, we assemble the voxel signal values into a 3D image matrix by Eq. (8). By repeating the multivoxel BOLD fMRI simulations for a range of T_{E} setting (T_{E} = 0:2:60 ms with an increment of 2 ms), for the different voxel sizes (128 × 128 × 128, 64 × 64 × 64, 32 × 32 × 32 micron^{3}), and for different vessel sizes (radii = 2 and 4 micron), we obtain a collection of MR magnitude images and phase images; from which we may observe the effect of MRI technology on the imaging performance of BOLD fMRI with respect to the echo time T_{E}, the image resolution, and the vessel size. In particular, we are concerned with the magnitudevssusceptibility correlation (corrA in Eq. (10)) and the phasevsfieldmap correlation (corrP in Eq. (11)). In what follows, we present the simulation results via figures.
Figure 5 shows a 3D geometry for a NAB embedded in a cortical FOV, in which the cortical vasculature is simulated with a mixture of randomly populated cylinders (radii = {2,4,8}micron). In order to show the effect of different vessel sizes, we carry out the simulations for monosized cortical vasculatures (with radii = 2 and 4 micron separately, under the condition of bfrac = 2%). Using monosized vasculature allows us to address the vesselsize effect on the BOLD fMRI signals (the BOLD fMRI nonlinearity due to large vessels [35]).
With 1micron digital grid resolution, we originally represent a FOV with large support matrix in size of 2048 × 2048 × 2048. After the fieldmap calculation, we partition the fieldmap into multivoxel image arrays for three voxel sizes: 16 × 16 × 16 matrix (voxel size: 128 × 128 × 128 micron^{3}), 32 × 32 × 32 matrix (voxel size: 64 × 64 × 64 micron^{3}), and 64 × 64 × 64 matrix (voxel size 32 × 32 × 32 micron^{3}). With these three voxel sizes, we demonstrate the effect of image resolution on BOLD fMRI.
2 μmradius vessels  4 μmradius vessels  

T_{E} = 1 ms  T_{E} = 30 ms  T_{E} = 1 ms  T_{E} = 30 ms  
Voxel size(μm^{3})  corrA  corrP  corrA  corrP  corrA  corrP  corrA  corrP 
32 × 32 × 32  0.867  1.000  0.885  0.998  0.834  1.000  0.854  0.998 
64 × 64 × 64  0.883  1.000  0.899  0.999  0.871  1.000  0.887  0.999 
128 × 128 × 128  0.892  1.000  0.907  0.999  0.898  1.000  0.910  0.999 
Discussion
In this paper, we propose a magnetic susceptibility expression formula for a linear neurovascular coupling model (accounting for the neurovascular aspects: CBF, CBV, CMRO_{2}, FOV, and NAB) and a computational model for volumetric BOLD fMRI simulation. Overall, our goal is to numerically examine the principles of MRIbased neuroimaging: implementing the forward imaging from a neuronal origin (a numerical NAB), to vascular response (a numerical susceptibility distribution Δχ, which also serves as the vascular origin of T2*MRI), to complex multivoxel MR image formation (MR magnitude images), and then rendering the backward mappings (as designated by "A ~ NAB", "A~Δχ" and "P~ΔB" in Figure 2). In this paper, we focus on the mapping between the MR magnitude image and its vascular origin. Meanwhile, we also show the mapping between the MR phase image and the fieldmap.
We decompose the overall BOLD fMRI model into two modules in Figure 2: "neurophysiology" and "MRI technology". Due to the complexity of biophysiology and neurology, the neurovascularcoupled BOLD process is not completely understood. In contrast, the MRI technology is well understood and developed. We propose a linear neurovascular coupling formula to express a BOLD state in terms of susceptibility perturbation in Eq. (4), which accounts for the main aspects of BOLD process as parameterized by CBV, CBF, CMRO_{2}, and NAB. It should be mentioned that the linear model neurovascular coupling model has the experimental supports [27, 28, 31, 32]. For numerical simulation, we need to express the dynamics of neurovascular coupling process and BOLD fMRI through the digitization of the involved parameters at discrete time points (snapshots). It is straightforward to extend our snapshot simulation method to dynamic BOLD fMRI simulation by repeating the snapshot imaging at a series of time points.
The cortical vasculature mainly consists of capillaries with radii in a range from 2 to 10 micron. Voxelization may suppress the intravoxel vasculature details (due to voxel size > > vessel size), consequently, we cannot discern the vascular structures in Figures 1, 6 through 8. Our simulation results show that, the magnitudevssusceptibility correlation decreases for larger vessel sizes, which may in part explain the experimentally observed BOLD nonlinearity due to large vessel effect [35]. In general, the magnitudevssusceptibility correlation takes on a correlation coefficient of 0.90 for an image resolution about 100micron. The discrepancy between the MR magnitude and the BOLD susceptibility map (corrA≠1) may be partially understood from the nonlinearity of BOLD fMRI [13, 35]. We also demonstrate the phasevsfieldmap correlation. The results show that, the corrP is close to the maximum 1.00 with respect to T_{E} (< 30 ms). However the corrP tends to drop for long T_{E} and high resolution (see Figure 9(b) and Table 2 and refer to [20] for more relevant results and explanations).
In our simulation implementation, the cortex FOV is represented by a large 3D support matrix in size of 2048 × 2048 × 2048 (with 1micron grid resolution) in order to depict small capillaries (digital geometry requirement). The output MR images are represented in much smaller multivoxel matrices (via voxelization): 16 × 16 × 16, 32 × 32 × 32, and 64 × 64 × 64 for three image resolutions. It is reminded that the fieldmap should be calculated from the finelygridded susceptibility map (in the large support matrix) as a whole by a FT technique (in Eq. (5)), rather than trying to divide the large matrix into smaller submatrices. The 3D FFT on a large matrix (that is too large to be processed as a whole in computer memory) is implemented by a scheme of 2D + 1D decomposition [36]. By performing the 3D FFT on the large matrix, we can account for longdistance magnetic influence of susceptibility magnetization [22, 23]. Indeed the multivoxel BOLD fMRI simulation demands a proliferated computation load in comparison with the single BOLD voxel signal simulation. Specifically, the volumetric simulation on an image matrix of 32 × 32 × 32 with voxel size of 64 × 64 × 64 (necessary for digital geometrical depiction of intravascular structure) requires a large matrix of 2048 × 2048 × 2048 (with 2048 = 64 × 32). Considering that the volumetric BOLD fMRI simulation may provide the spatial features among voxels that is not available from individual voxel signals, our simulation effort is worthwhile.
The simulation results for magnitudebased and phasebased backward mappings are graphically shown in Figure 9, from which we can observe the following phenomenon: 1) The magnitudevssusceptibility correlation (corrA) increases with respect to T_{E} (shown in a range from 0 to 60 ms), whereas the phasevsfieldmap correlation (corrP) decreases for T_{E} > 30 ms; 2) Large voxel size (or low image resolution) cause both corrA and corrP decrease (noticeable corrP decrease for T_{E} > 30 ms); 3) Large vessels cause both corrA and corrP drop.
In addition to the study on the backward mapping between the MR magnitude image and its vascular origin, we also provide the backward mapping between the MR phase image and the fieldmap in a measure of corrP. Our simulation results show that corrP ≈ 1.00, which indicates that the MR phase image conforms very well with the fieldmap. This suggests the possibility of susceptibility reconstruction by computed inverse MRI [13, 26], which is an active ongoing research topic.
Numerical simulation provides a powerful tool to look into a digital imaging system as long as all the involved parameters can be numerically characterized. In light of numerical method, a parameter must be digitized for computation. The parameter digitization does not necessarily require an analytic formula, it may be fulfilled in any way. For example, the BOLD susceptibility perturbation over a FOV in Figure 4 can be assigned a NABmodulated spatial distribution at a snapshot time, and the dynamics can be numerically characterized by an empirical timecourse in a lookup table with respect to a discrete time (the voxel timecourses in Figure 4 are not necessarily analytically describable). Therefore, our computational model can be extended to accommodate a nonlinear neurovascularcoupled BOLD fMRI process over a broad physiological range provided that all the involved parameters can be numerically determined.
Conclusions
In this report we propose a computational model for numerically simulating volumetric BOLD MRI with a magnetic susceptibility expression for a linear neurovascular coupling state. The forward procedure for this model includes: 1) defining a neuronal origin (a numerical NAB) in a cortex region(FOV); 2) expressing the NABmodulated vascular response by a spatial distribution of intravascular blood biomagnetic susceptibility perturbation; and 3) calculating the susceptibilityinduced fieldmap by accounting for the magnetization of all vascular blood in the FOV; 4) calculating the multivoxel MR image by intravoxel dephasing integration. Upon the completion of the forward simulation, we compare the output MR magnitude image with the predefined susceptibility map in a numerical measure of magnitudevssusceptibility correlation, thereby quantitatively examining the reproducibility of the vascular origin by the MR magnitude image. Meanwhile, we can also compare the output MR phase image with the precomputed fieldmap in terms of phasevsfieldmap correlation, showing that the MR phase image conforms with the fieldmap very well.
Based on our simulation results, we conclude that, 1) the magnitude image of BOLD fMRI can approximately, but not exactly, represent the vascular origin; and 2) the phase image conforms very well with the fieldmap. Considering that the vascular response to the neuronal origin is subject to a neurovascular coupling process, we can infer that the mapping between the MR magnitude image and neuronal origin suffers more mismatch than that between the MR magnitude image and the vascular origin, even if under a linear neurovascular coupling model (because the vasculature exerts a spatial modulation that imposes additional mismatch). The volumetric computational model provides a general framework for simulating many neurovascular, physiological, and biophysical aspects. It can be extended to accommodate a nonlinear neurovascular coupling process over a broad physiological range if the magnetic susceptibility expression is numerically available (not necessarily formulable). Our simulation results pose a caveat to the MRIbased neuroimaging and brain mapping study: the MR magnitude image is not an exact reproduction which may in part explain the nonlinearity of BOLD fMRI. In future research we plan to to look into the overall nonlinearity of BOLD fMRI by incorporating the intrinsic nonlinear neurovascular coupling process.
Abbreviations
 BOLD:

Blood oxygenation level dependent
 MR:

Magnetic resonance
 MRI:

Magnetic resonance imaging
 fMRI:

Functional magnetic resonance imaging
 FOV:

Field of view
 NAB:

Neuroactive blob
 CBV:

Cerebral blood volume
 CBF:

Cerebral blood flow
 CMRO_{2}:

Cerebral metabolic rate of oxygen
 bfrac :

Blood volume fraction.
Declarations
Acknowledgements
This research was supported by the NIH (1R01EB006841, 1R01EB005846), NSF (0612076), and the MRN internal grant 6003154.
Authors’ Affiliations
References
 Boxerman JL, Bandettini PA, Kwong KK, Baker JR, Davis TL, Rosen BR, Weisskoff RM: The intravascular contribution to fMRI signal change: Monte Carlo modeling and diffusionweighted studies in vivo. Magn Reson Med. 1995, 34 (1): 410. 10.1002/mrm.1910340103.View ArticlePubMed
 Boxerman JL, Hamberg LM, Rosen BR, Weisskoff RM: MR contrast due to intravascular magnetic susceptibility perturbations. Magn Reson Med. 1995, 34 (4): 555566. 10.1002/mrm.1910340412.View ArticlePubMed
 Ogawa S, Menon RS, Tank DW, Kim SG, Merkle H, Ellermann JM, Ugurbil K: Functional brain mapping by blood oxygenation leveldependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model. Biophys J. 1993, 64 (3): 803812. 10.1016/S00063495(93)814413.PubMed CentralView ArticlePubMed
 Weisskoff RM, Zuo CS, Boxerman JL, Rosen BR: Microscopic susceptibility variation and transverse relaxation: theory and experiment. Magn Reson Med. 1994, 31 (6): 601610. 10.1002/mrm.1910310605.View ArticlePubMed
 Yablonskiy DA, Haacke EM: Theory of NMR signal behavior in magnetically inhomogeneous tissues: the static dephasing regime. Magn Reson Med. 1994, 32 (6): 749763. 10.1002/mrm.1910320610.View ArticlePubMed
 Pathak AP, Ward BD, Schmainda KM: A novel technique for modeling susceptibilitybased contrast mechanisms for arbitrary microvascular geometries: the finite perturber method. Neuroimage. 2008, 40 (3): 11301143. 10.1016/j.neuroimage.2008.01.022.PubMed CentralView ArticlePubMed
 Arthurs OJ, Boniface S: How well do we understand the neural origins of the fMRI BOLD signal?. Trends Neurosci. 2002, 25 (1): 2731. 10.1016/S01662236(00)019950.View ArticlePubMed
 Attwell D, Iadecola C: The neural basis of functional brain imaging signals. Trends Neurosci. 2002, 25 (12): 621625. 10.1016/S01662236(02)022646.View ArticlePubMed
 Buxton RB: Interpreting oxygenationbased neuroimaging signals: the importance and the challenge of understanding brain oxygen metabolism. Front Neuroenergetics. 2010, 2: 8PubMed CentralPubMed
 Buxton RB, Frank LR: A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation. J Cereb Blood Flow Metab. 1997, 17 (1): 6472.View ArticlePubMed
 Zhang N, Liu Z, He B, Chen W: Noninvasive study of neurovascular coupling during graded neuronal suppression. Journal of cerebral blood flow and metabolism: official journal of the International Society of Cerebral Blood Flow and Metabolism. 2008, 28 (2): 280290.View Article
 Chavhan GB, Babyn PS, Thomas B, Shroff MM, Haacke EM: Principles, techniques, and applications of T2*based MR imaging and its special applications. Radiographics: a review publication of the Radiological Society of North America, Inc. 2009, 29 (5): 14331449.View Article
 Chen Z, Calhoun VD: Two pitfalls of BOLD fMRI magnitudebased neuroimage analysis: nonnegativity and edge effect. J Neurosci Methods. 2011, 199 (2): 363369. 10.1016/j.jneumeth.2011.05.018.PubMed CentralView ArticlePubMed
 Howseman AM, Bowtell RW: Functional magnetic resonance imaging: imaging techniques and contrast mechanisms. Philos Trans R Soc Lond B Biol Sci. 1999, 354 (1387): 11791194. 10.1098/rstb.1999.0473.PubMed CentralView ArticlePubMed
 Moon CH, Fukuda M, Park SH, Kim SG: Neural interpretation of blood oxygenation leveldependent fMRI maps at submillimeter columnar resolution. J Neurosci. 2007, 27 (26): 68926902. 10.1523/JNEUROSCI.044507.2007.View ArticlePubMed
 Marques JP, Bowtell RW: Using forward calculations of the magnetic field perturbation due to a realistic vascular model to explore the BOLD effect. NMR Biomed. 2008, 21 (6): 553565. 10.1002/nbm.1224.View ArticlePubMed
 Martindale J, Kennerley AJ, Johnston D, Zheng Y, Mayhew JE: Theory and generalization of Monte Carlo models of the BOLD signal source. Magn Reson Med. 2008, 59 (3): 607618. 10.1002/mrm.21512.View ArticlePubMed
 Menon RS: Simulation of BOLD phase and magnitude changes in a voxel. Proc Intl Soc Mag Reson Med. 2003, 11: 1719
 Chen Z, Calhoun V: A computational multiresolution BOLD fMRI model. IEEE Trans BioMed Eng. 2011, 58 (10): 29952999.PubMed CentralView ArticlePubMed
 Chen Z, Calhoun VD: Magnitude and Phase Behavior of Multiresolution BOLD signal. Concepts Magn Reson. 2010, 37B (3): 129135. 10.1002/cmr.b.20164.View Article
 Chen Z, Chen Z, Calhoun VD: Multiresolution voxel decomposition of complexvalued BOLD signals reveals phasor turbulence. Orlando: In: SPIE Medical Imaging;:111. 16 March 79613y; Orlando, FL: SPIE; 79613y
 Chen Z, Chen Z, Calhoun VD: Voxel magnetic field disturbance from remote vasculature in BOLD. fMRI.In:. Orlando: SPIE Medical Imaging; 112. 16 March 79613x; Orlando, FL: SPIE; 79613x
 Chen Z, Caprihan A, Calhoun VD: Effect of surrounding vasculature on intravoxel BOLD signal. Med Phys. 2010, 37 (4): 17781787. 10.1118/1.3366251.PubMed CentralView ArticlePubMed
 Koch KM, Papademetris X, Rothman DL, de Graaf RA: Rapid calculations of susceptibilityinduced magnetostatic field perturbations for in vivo magnetic resonance. Phys Med Biol. 2006, 51 (24): 63816402. 10.1088/00319155/51/24/007.View ArticlePubMed
 Salomir R, de Senneville BD, Moonen CTW: A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts Magn Reson. 2003, B19: 2634.View Article
 Chen Z, Calhoun V: Computed inverse MRI for magnetic susceptibility map reconstruction. J Comp Assit Tomography. 2012, 36 (2): 265274. 10.1097/RCT.0b013e3182455cab. March/April 2012View Article
 Arthurs OJ, Williams EJ, Carpenter TA, Pickard JD, Boniface SJ: Linear coupling between functional magnetic resonance imaging and evoked potential amplitude in human somatosensory cortex. Neuroscience. 2000, 101 (4): 803806. 10.1016/S03064522(00)00511X.View ArticlePubMed
 Hoge RD, Atkinson J, Gill B, Crelier GR, Marrett S, Pike GB: Linear coupling between cerebral blood flow and oxygen consumption in activated human cortex. Proc Natl Acad Sci USA. 1999, 96 (16): 94039408. 10.1073/pnas.96.16.9403.PubMed CentralView ArticlePubMed
 Logothetis NK: The neural basis of the bloodoxygenleveldependent functional magnetic resonance imaging signal. Philos Trans R Soc Lond B Biol Sci. 2002, 357 (1424): 10031037. 10.1098/rstb.2002.1114.PubMed CentralView ArticlePubMed
 Ugurbil K, Toth L, Kim DS: How accurate is magnetic resonance imaging of brain function?. Trends Neurosci. 2003, 26 (2): 108114. 10.1016/S01662236(02)000395.View ArticlePubMed
 Zhang N, Yacoub E, Zhu XH, Ugurbil K, Chen W: Linearity of bloodoxygenationlevel dependent signal at microvasculature. Neuroimage. 2009, 48 (2): 313318. 10.1016/j.neuroimage.2009.06.071.PubMed CentralView ArticlePubMed
 Norris DG: Principles of magnetic resonance assessment of brain function. Journal of magnetic resonance imaging: JMRI. 2006, 23 (6): 794807. 10.1002/jmri.20587.View ArticlePubMed
 Uludag K, MullerBierl B, Ugurbil K: An integrative model for neuronal activityinduced signal changes for gradient and spin echo functional imaging. Neuroimage. 2009, 48 (1): 150165. 10.1016/j.neuroimage.2009.05.051.View ArticlePubMed
 Marques JP, Bowtell R: Application of a fourierbased method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts Magn Reson. 2005, 6578. B25:
 Zhang N, Zhu XH, Chen W: Investigating the source of BOLD nonlinearity in human visual cortex in response to paired visual stimuli. Neuroimage. 2008, 43 (2): 204212. 10.1016/j.neuroimage.2008.06.033.PubMed CentralView ArticlePubMed
 Chen Z, Ning R: Breast volume denoising and noise characterization by 3D wavelet transform. Comput Med Imaging Graph. 2004, 28 (5): 235246. 10.1016/j.compmedimag.2004.04.004.View ArticlePubMed
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14712342/12/8/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.