We motivate our approach by looking at a typical neuroimage as shown in Figure , 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 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 ; 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.
Figure 1 Illustration of a local neuroactive blob (NAB) observed at cerebral cortex by fMRI experiment. Conventionally, a NAB is inferred from an fMRI data and is considered as a local neuronal activity. However this convention has not been numerically verified (more ...)
The overall diagram for the computational model of neurovascular coupling and BOLD fMRI is shown in Figure , 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 , 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 ). 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 .
Figure 2 Overall diagram of the computational model for neurovascular coupling and BOLD fMRI. It is decomposed into a neurophysiology module (upper dashed box) and a MRI technology module (down dashed box). The linear neurovascular coupling model accounts for (more ...)
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 D0 × D0 × D0) by
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), (x0, y0, z0) 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 can be formulated by a combination of several regular Gaussian-shaped NAB primitives that are different in (x0, y0, z0) 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 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).
Figure 3 A 2D illustration of a localized neurovascular coupling model. The local neuroactivity is represented by a ball-shaped neuroactive blob (NAB), which confines the locality and the graded strength in a cortical FOV. The red and blue lines illustrate the (more ...)
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,
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 (CMRO2). 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
= 0.27 × 4π
ppm (SI unit) represents the magnetic susceptibility change between deoxyhemoglobin and oxyhemoglobin, which has been used for BOLD signal simulation [17
]. 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
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, CMRO2 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 ).
Parameters and settings for numerical simulations
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 CMRO2. 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 (Δχ
(1-Y)) that has been accepted for single voxel BOLD signal simulation [1
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, CMRO2, and NAB on the BOLD susceptibility perturbation (which is detected by T2*MRI). In Figure , 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 t0 through t5 in Figure ), 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.
Figure 4 Illustration of dynamic magnetic susceptibility perturbation in a NAB-weighted cortical region. It is shown that the voxel at the central NAB region produces the maximum susceptibility perturbation (in comparison with the voxels outside the central NAB), (more ...)
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 B0
), called the fieldmap henceforth, by [24
where FT and IFT stand for Fourier transform and inverse Fourier transform respectively, (kx, ky, kz) 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 shows a 3D FOV of size D0
, which is filled with random vascular networks (the cortical vasculature generation technique has been reported previously [20
]). The FOV is partitioned into voxels (voxel size = d0
), 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
Figure 5 A typical geometry of cortical FOV that is filled with vasculature and encloses a local neuroactivity blob (NAB in red). The fMRI-detectable neurovascular coupling state is expressed as a NAB-weighted intravascular blood magnetic susceptibility perturbation (more ...)
where * denotes the convolution, and V [xn, yn, zn] 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[xn, yn, zn], Δχ[xn, yn, zn], and ΔB[xn, yn, zn], 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, TE
) = γ
(x, y, z
, where γ
is the gyromagnetic ratio. It is noted that the phase angle is different from the field value by a constant factor γ
(Larmor law). Due to the finite dimension of a voxel (for example, d0
= 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 d0
voxel at discrete position [xn, yn, zn
], its BOLD signal is formed by
where the echo time TE is explicitly retained to remind of the TE 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 TE), we assemble the voxel signal values into a 3D MR matrix according to the voxelization scheme in Eq. (6) by
which is complex-valued and explicitly TE-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 TE. The magnitude loss map and phase accumulation map for a given TE are calculated by
= 0]|denotes the non-decay initial magnitude image and
= 0] the initial phase image. Henceforth, we will denote the MR magnitude and phase images as A[xn
] and P[xn
In the forward BOLD fMRI calculation, we obtain a pair of MR magnitude image and phase image by Eq. (9). Considering the NAB[xn, yn, zn] (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 . 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 .
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[xn, yn, zn] and the predefined BOLD susceptibility map Δχ[xn, yn, zn] by a spatial correlation coefficient as defined by
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[xn
]→0 for TE
→ 0), a meaningful corrA(TE
) should be evaluated at a relative long TE
> 0). However, a long TE
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 , which connotes the magnitude-vs-susceptibility correlation.
Likewise, we can measure the spatial correlation between the MR phase image P[xn, yn, zn; TE] and the fieldmap ΔB[xn, yn, zn] by
Likely, corrP is a numerical measure of the phase-vs-fieldmap correlation, a backward mapping designated by "P~ΔB" in Figure .
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.