The blood oxygenation-level 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 neurovascular-coupled 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 susceptibility-induced 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 susceptibility-induced fieldmap (spatial correlation ≈ 1.00).

Blood oxygenation-level 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 neuroactivity-induced 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 susceptibility-expressed 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 neurovascular-coupled 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 multi-voxel 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 susceptibility-expressed 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 intra-FOV 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 complex-valued 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 susceptibility-expressed 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

We motivate our approach by looking at a typical neuroimage as shown in Figure 1, which can be modeled as a population of local neuroactivity blobs (NAB). A spatially localized NAB is inferred from an fMRI dataset and is interpreted as (or related to) the underlying neuronal activity (neuronal origin). Because this inference is based on MR magnitude images, it is important to evaluate how well a MR magnitude image can represent a snapshot of the neuronal origin. Since the T2*MRI detection is a typical noninvasive 3D imaging modality, it is possible (virtually necessary) to quantitatively assess the neuronal origin inference from MR images by following the medical imaging principle: predefining an input source, predicting its output image by numerical simulation (or phantom experimental verification when applicable), and then comparing the obtained output image with the predefined input source. Accordingly, we address the NAB inference issue in Figure 1 by numerically simulating the neurovascular coupling process and the BOLD fMRI detection technology by the following steps: 1) delineating a cortical FOV that is ample enough to enclose a NAB, as demonstrated in Figure 1; 2) defining a susceptibility-expressed BOLD state under the NAB-modulated neurovascular coupling; and 3) carrying out the multivoxel BOLD fMRI simulation for a predefined susceptibility source. Upon the completion of BOLD fMRI simulation, we compare the MR magnitude image with the predefined susceptibility source in a measure of spatial correlation. A key benefit of our approach is we produce a full 3D image in form of a multivoxel matrix, in contrast to the many previous studies which focused only on one or a few voxels.

The overall diagram for the computational model of neurovascular coupling and BOLD fMRI is shown in Figure 2, which consists of two cascaded modules: 1) the "neurophysiology" module for the neurovascular coupling process from a neuronal origin to its vascular response, which should be phenotypically expressed in terms of biomagnetic susceptibility distribution for the purpose of MRI detection, and 2) the "MRI technology" module for imaging a susceptibility-expressed BOLD state by T2*MRI. According to this two-module BOLD fMRI model, the neuronal origin is a NAB-expressed neuroactivity (in response to an external stimulus, not shown), and the vascular response is a result of neurovascular coupling (numerically expressed in form of Δχ), which in turn serves as the vascular origin of T2*MRI detection. Conventionally, the goal of brain mapping and neuroimaging is to depict the neuroactivity origin from MR magnitude images, as diagrammed by a grey double-directed arrow for the backward mapping as designated by "A ~ NAB". Based on the two-module model in Figure 2, the overall backward mapping of "A ~ NAB" can be decomposed into two steps: the first step is from the MR magnitude image to the vascular origin, as denoted by a black double-directed arrow and designated by "A~Δχ", the second step is from the vascular origin to the neuronal origin (not shown in Figure 2). In this paper, we will look into the first step backward mapping of "A~Δχ" by numerical T2*MRI simulation, and simplify the second step backward mapping by a linear neurovascular coupling model. Meanwhile, we also show another backward mapping: from a MR phase image to fieldmap as diagrammed by a black double-directed arrow as designated by "P~ΔB" in Figure 2.

Neurovascular coupling formulation

It is known in neurobiology and neurophysiology that a neuroactivity is accompanied by a complicated process of cellular, metabolic, and vascular processes. For simplification of computational implementation, we express the spatial distribution for the functional parcellation of a neuroactivity by a 3D Gaussian-shaped NAB embedded in a cortical FOV (with size D_{0} × D_{0} × D_{0}) by

(1)

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 non-spherical blob like the one in Figure 1 can be formulated by a combination of several regular Gaussian-shaped 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).

In response to the neuroactivity in a NAB, the neurovascular unit regulates the cerebral microcirculation by vasodilation (increase in CBV) and a burst of blood inflow (increase in CBF). The intravascular blood oxygenation level is subject to change due to neuromodulation in the NAB, primarily occurring in the capillaries and venules. Figure 3 illustrates a neurovascular coupled blood oxygen delivery in capillary bed confined in a NAB. It is shown that the BOLD activity is spatially localized in a cortex region (NAB confinement) where the central region undergoes the maximal activity, and there is no neuronal activity outside NAB. The blood biomagnetic susceptibility perturbation is determined by the spatial modulation between the NAB distribution and the vasculature geometry. The intravascular blood volume (CBV) is determined by the vasculature geometry at a snapshot pose, and the dynamic intravascular blood flow (CBF) is determined by the intracranial blood pressure, vasculature geometry and blood fluid mechanics. All in all, for numerical simulation purpose, all these neurovascular parameters should be numerically expressed in terms of magnetic susceptibility responses (detectable by MRI technology).

At a snapshot time, the cortical vasculature geometry determines the intravascular blood volume, that only takes up a fractional space of the cortical FOV, as described by blood volume fraction (bfrac). For human cerebral cortex, bfrac ≈ 2% [33]. For numerical simulation of 3D vasculature-laden FOV, we express the vasculature configuration by a binary volume that is randomly generated under the control of bfrac. That is,

(2)

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 (1-Y(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.

Based on the neurovascular-coupled blood biomagnetic perturbation mechanism, we propose a biomagnetic susceptibility expression formula for a linear neurovascular coupling model by

(3)

where χ_{do} = χ_{deoxy}-χ_{oxy} = 0.27 × 4πppm (SI unit) represents the magnetic susceptibility change between deoxyhemoglobin and oxyhemoglobin, which has been used for BOLD signal simulation [17, 23]. The total susceptibility χ_{total} includes contributions from both intravascular blood and extravascular tissue parenchyma. The susceptibility distribution of a selected baseline state is denoted by χ_{base}. In reference to χ_{base}, we can characterize a BOLD susceptibility perturbation state by Δχ in Eq. (3). Suppose the BOLD fMRI system is a linear digital imaging system, then the behavior of MR image change in reference to its baseline state can represent the intravascular susceptibility perturbation Δχ(x, y, z, t). In other words, the linear BOLD fMRI model allows us to infer the intravascular BOLD susceptibility perturbation by observing the corresponding MR image change (in reference to their respective baselines). From Eq. (3), we can express the NAB-modulated BOLD susceptibility perturbation state by

(4)

This is a computational model for linear neurovascular coupling, which provides a mathematical formula for numerically expressing the NAB-modulated BOLD response process in a vasculature-filled FOV. Specifically, CMRO_{2} is accounted for by Y(t), CBF and CBV by bfrac(t) that is implicitly embodied in V(x, y, z, t) via vasculature configuration (during vasculature geometry generation under the condition of bfrac(t) in Eq. (2)), and the local neuroactivity by NAB(x, y, z) which confines the neuroactivity extent and defines a graded neuroactivity strength. Usually, the blood magnetism parameter χ_{do} and the blood physiology parameter Hct assume for normal blood, which are experimentally determined constants (see Table 1).

For the neurovascular coupling formula in Eq. (4), we need to point out following aspects:

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 (1-Y(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 1-Y = 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⋅(1-Y)) that has been accepted for single voxel BOLD signal simulation [1, 3, 17, 23].

Overall, a multivoxel BOLD fMRI simulation requires a predefined magnetic susceptibility distribution as the input source of T2*MRI. The susceptibility expression for a neurovascular coupling process plays a bridge between the neuroscience and MRI technology. Although the neurovascular coupling process is not fully understood so far, we propose a linear spatiotemporal modulation model (in Eq. (4)) that allows us to look into the effects of CBF, CBV, CMRO_{2}, and NAB on the BOLD susceptibility perturbation (which is detected by T2*MRI). In Figure 4, we illustrate a dynamic susceptibility perturbation by susceptibility timecourses of 3 voxels in a NAB, which shows that the susceptibility perturbation strength is spatially weighted by the NAB (maximum at the center and reduced toward boundary). According to Eq. (4), the numerical characterization of neurovascular-coupled 4D dynamic susceptibility perturbation Δχ(x, y, z, t) involves the following conditions: defining a local neuroactivity (a numerical NAB that is not necessarily in an analytic formula), generating a cortical vasculature (V(x, y, z, t)) under the control of bfrac(t), and assigning a value to the blood oxygenation level Y(t) (= [0, 1]). It is mentioned that if the digitization of Δχ(x, y, z, t) are experimentally or empirically available for a train of discrete time points (such as at the ticks t_{0} through t_{5} in Figure 4), it is not necessary to seek the analytic formula in Eq. (4). In the following, we will implement the numerical T2*MRI simulation for acquiring a complex MR image from a given susceptibility distribution at a snapshot time, denoted by Δχ(x, y, z), thus omitting the time variable t.

Forward BOLD fMRI simulation

From magnetism perturbation to fieldmap establishment

Given a snapshot of neurovascular-coupled BOLD state, Δχ(x, y, z), we can calculate its magnetization field distribution (resulting from the blood magnetization in a main field B_{0}), called the fieldmap henceforth, by [24, 25, 34]

(5)

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

We simulate the cortical FOV by filling it with vasculature. Figure 5 shows a 3D FOV of size D_{0} × D_{0} × D_{0}, which is filled with random vascular networks (the cortical vasculature generation technique has been reported previously [20, 23]). The FOV is partitioned into voxels (voxel size = d_{0} × d_{0} × d_{0}), thereby we can represent the FOV in a small array of voxels, in which each voxel can be assigned a value by intravoxel average. The process of spatial partition into voxels and intravoxel average is called voxelization. For the vasculature V(x, y, z) in the FOV, the voxelization is expressed by

(6)

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

Exposed to an inhomogeneous fieldmap, a proton precesses with a phase angle Δϕ(x, y, z, T_{
E
}) = γ · ΔB(x, y, z)·T_{
E
}, where γ is the gyromagnetic ratio. It is noted that the phase angle is different from the field value by a constant factor γ⋅T_{E} (Larmor law). Due to the finite dimension of a voxel (for example, d_{0} = 128 micron in our simulation), its voxel signal is formed by a vector sum of all spin packets (or isochromats) inside the voxel, called intravoxel dephasing integration [19]. For a d_{0} × d_{0} × d_{0} voxel at discrete position [x_{
n
}, y_{
n
}, z_{
n
}], its BOLD signal is formed by

(7)

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

After calculating the voxel signals for all voxels in the FOV (with a specific T_{E}), we assemble the voxel signal values into a 3D MR matrix according to the voxelization scheme in Eq. (6) by

(8)

which is complex-valued and explicitly T_{E}-dependent. Since the MR magnitude image contrast is due to a spatial distribution of voxel signal decay, we are concerned with the magnitude loss. For the phase image, we are concerned with the phase angle accumulation during the period of T_{E}. The magnitude loss map and phase accumulation map for a given T_{E} are calculated by

(9)

Where |C[x_{n}, y_{n}, z_{n}; T_{E} = 0]|denotes the non-decay 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.

It is known that the underlying source of BOLD fMRI is the susceptibility-expressed BOLD state, which is highly dependent upon the vasculature configuration in the FOV. To reduce the effect of vessel randomness on our numerical simulation, we propose to measure the similarity of the MR magnitude image A[x_{n}, y_{n}, z_{n}] and the predefined BOLD susceptibility map Δχ[x_{n}, y_{n}, z_{n}] by a spatial correlation coefficient as defined by

(10)

Where cov(x, y) denotes the covariance between vector x and y, std(x) the standard deviation, and ":" a nD-to-1D operator (as used by Matlab language) that reorders a high-dimensional 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 magnitude-vs-susceptibility correlation.

Likewise, we can measure the spatial correlation between the MR phase image P[x_{n}, y_{n}, z_{n}; T_{E}] and the fieldmap ΔB[x_{n}, y_{n}, z_{n}] by

(11)

Likely, corrP is a numerical measure of the phase-vs-fieldmap correlation, a backward mapping designated by "P~ΔB" in Figure 2.

In summary, our computational BOLD fMRI model can be used for backward mappings: magnitude-vs-susceptibility correlation (corrA) and phase-vs-fieldmap 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 neurovascular-coupled 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 Gaussian-shaped NAB, a vasculature-laden 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 magnitude-vs-susceptibility correlation (corrA in Eq. (10)) and the phase-vs-fieldmap 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 vessel-size effect on the BOLD fMRI signals (the BOLD fMRI nonlinearity due to large vessels [35]).

With 1-micron 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.

Figure 6 shows the NAB-profiled susceptibility perturbation map (calculated by Eq. (4)) and the subsequent fieldmap (calculated by Eq. (5)). Figure 7 shows the MR magnitude and phase images (displayed with z-slices from FOV surface to the center surface) calculated at T_{E} = 30 ms for voxel size = 64 × 64 × 64 micron^{3}. It is noted that the magnitude image resembles the susceptibility map to a great extent, and the phase image replicates the fieldmap very well. Figure 8 shows the T_{E} dependence of the MR magnitude and phase images (with the central z-slice).

With the MR magnitude and phase images generated for a variety of parameter settings, we calculate the magnitude-vs-susceptibility correlation and the phase-vs-fieldmap correlation by Eqs. (10) and (11), respectively. The results are shown in Figure 9 for T_{E} in a range of 0 to 60 ms and three different image resolutions (see the legend). From Figure 9 we can see that the magnitude-vs-susceptibility correlation increases with respect to T_{E} and voxel size; on contrary, the phase-vs-fieldmap correlation decreases with respect to T_{E} and voxel size.

In Table 2 we present the correlation coefficients for two vasculatures (radii = 2 and 4 micron), three image resolutions (voxel size = 128 × 128 × 128, 64 × 64 × 64, and 32 × 32 × 32 micron^{3}), and two selected echo times (T_{E} = 1 ms and 30 ms).

Table 2

Spatial matching results in terms of spatial correlation coefficients corrA and corrP calculated by Eqs.(10) and (11) with the setting of B_{0} = 3 T, bfrac = 2%, radii = 2 and 4 μm, T_{E} = 1 and 30 ms for three image resolutions (voxel sizes)

2- μm-radius vessels

4- μm-radius 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 MRI-based 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 neurovascular-coupled 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 magnitude-vs-susceptibility 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 magnitude-vs-susceptibility correlation takes on a correlation coefficient of 0.90 for an image resolution about 100-micron. 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 phase-vs-fieldmap 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 1-micron 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 finely-gridded 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 long-distance 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 magnitude-based and phase-based backward mappings are graphically shown in Figure 9, from which we can observe the following phenomenon: 1) The magnitude-vs-susceptibility correlation (corrA) increases with respect to T_{E} (shown in a range from 0 to 60 ms), whereas the phase-vs-fieldmap 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 NAB-modulated 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 neurovascular-coupled 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 NAB-modulated vascular response by a spatial distribution of intravascular blood biomagnetic susceptibility perturbation; and 3) calculating the susceptibility-induced 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 magnitude-vs-susceptibility 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 phase-vs-fieldmap 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 MRI-based 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 6003-154.

Authors’ Affiliations

(1)

The Mind Research Network and LBERI

(2)

Department of Electrical and Computer Engineering

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 diffusion-weighted studies in vivo.Magn Reson Med 1995,34(1):4–10.PubMedView Article

Boxerman JL, Hamberg LM, Rosen BR, Weisskoff RM: MR contrast due to intravascular magnetic susceptibility perturbations.Magn Reson Med 1995,34(4):555–566.PubMedView Article

Ogawa S, Menon RS, Tank DW, Kim SG, Merkle H, Ellermann JM, Ugurbil K: Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model.Biophys J 1993,64(3):803–812.PubMedView Article

Weisskoff RM, Zuo CS, Boxerman JL, Rosen BR: Microscopic susceptibility variation and transverse relaxation: theory and experiment.Magn Reson Med 1994,31(6):601–610.PubMedView Article

Yablonskiy DA, Haacke EM: Theory of NMR signal behavior in magnetically inhomogeneous tissues: the static dephasing regime.Magn Reson Med 1994,32(6):749–763.PubMedView Article

Pathak AP, Ward BD, Schmainda KM: A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: the finite perturber method.Neuroimage 2008,40(3):1130–1143.PubMedView Article

Arthurs OJ, Boniface S: How well do we understand the neural origins of the fMRI BOLD signal?Trends Neurosci 2002,25(1):27–31.PubMedView Article

Attwell D, Iadecola C: The neural basis of functional brain imaging signals.Trends Neurosci 2002,25(12):621–625.PubMedView Article

Buxton RB: Interpreting oxygenation-based neuroimaging signals: the importance and the challenge of understanding brain oxygen metabolism.Front Neuroenergetics 2010, 2:8.PubMed

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):64–72.PubMedView Article

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):280–290.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):1433–1449.

Chen Z, Calhoun VD: Two pitfalls of BOLD fMRI magnitude-based neuroimage analysis: non-negativity and edge effect.J Neurosci Methods 2011,199(2):363–369.PubMedView Article

Howseman AM, Bowtell RW: Functional magnetic resonance imaging: imaging techniques and contrast mechanisms.Philos Trans R Soc Lond B Biol Sci 1999,354(1387):1179–1194.PubMedView Article

Moon CH, Fukuda M, Park SH, Kim SG: Neural interpretation of blood oxygenation level-dependent fMRI maps at submillimeter columnar resolution.J Neurosci 2007,27(26):6892–6902.PubMedView Article

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):553–565.PubMedView Article

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):607–618.PubMedView Article

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):2995–2999.PubMedView Article

Chen Z, Calhoun VD: Magnitude and Phase Behavior of Multiresolution BOLD signal.Concepts Magn Reson 2010,37B(3):129–135.View Article

Chen Z, Chen Z, Calhoun VD: Multiresolution voxel decomposition of complex-valued BOLD signals reveals phasor turbulence. Orlando: In: SPIE Medical Imaging;:1–11. 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; 1–12. 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):1778–1787.PubMedView Article

Koch KM, Papademetris X, Rothman DL, de Graaf RA: Rapid calculations of susceptibility-induced magnetostatic field perturbations for in vivo magnetic resonance.Phys Med Biol 2006,51(24):6381–6402.PubMedView Article

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:26–34.View Article

Chen Z, Calhoun V: Computed inverse MRI for magnetic susceptibility map reconstruction.J Comp Assit Tomography 2012,36(2):265–274. 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):803–806.PubMedView Article

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):9403–9408.PubMedView Article

Logothetis NK: The neural basis of the blood-oxygen-level-dependent functional magnetic resonance imaging signal.Philos Trans R Soc Lond B Biol Sci 2002,357(1424):1003–1037.PubMedView Article

Ugurbil K, Toth L, Kim DS: How accurate is magnetic resonance imaging of brain function?Trends Neurosci 2003,26(2):108–114.PubMedView Article

Zhang N, Yacoub E, Zhu XH, Ugurbil K, Chen W: Linearity of blood-oxygenation-level dependent signal at microvasculature.Neuroimage 2009,48(2):313–318.PubMedView Article

Norris DG: Principles of magnetic resonance assessment of brain function.Journal of magnetic resonance imaging: JMRI 2006,23(6):794–807.PubMedView Article

Uludag K, Muller-Bierl B, Ugurbil K: An integrative model for neuronal activity-induced signal changes for gradient and spin echo functional imaging.Neuroimage 2009,48(1):150–165.PubMedView Article

Marques JP, Bowtell R: Application of a fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility.Concepts Magn Reson 2005, (B25:):65–78.

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):204–212.PubMedView Article

Chen Z, Ning R: Breast volume denoising and noise characterization by 3D wavelet transform.Comput Med Imaging Graph 2004,28(5):235–246.PubMedView Article