**4. Discussion**

The data acquisition of BOLD fMRI is not analytically tractable due to the involvement of diversified parameters. The BOLD fMRI simulations provide a means to observe the effect of MRI transformations on the MR data acquisition; spatial distortions between the underlying magnetic source and MR images; and reproducibility of functional activity extraction from a 4D BOLD fMRI dataset.

Since MRI is designed to measure a magnetic field distribution, the BOLD fMRI only measures the χ-expressed BOLD response during a functional activity, a phenotypic numeric expression of a biophysiological brain functional activity in terms of tissue magnetism. It is believed that a functional activity causes IV blood magnetism change in terms of oxygenated and deoxy‐ genated blood magnetic susceptibility change. Therefore, our simulation begins with a configuration of magnetic source by a vasculature-laden VOI with a 3D χ source distribution. Through a spatiotemporal modulation by a predefined local neuroactive blob (numerically NAB(**r**)) and a task paradigm (numerically task(*t*)), we define a dynamic χ source to represent a χ-expressed BOLD activity (in Eq. (1)). It is pointed out in Eq. (1) that the BOLD χ response may incorporate the factors of cerebral metabolic rate of oxygen consumption (CMRO2), cerebral blood volume (CBV), and cerebral blood flow (CBF) through the parameters of *Y*(**t**) and *V*(**r**,*t*), thereby enabling the numeric simulations of MRI-detected BOLD activity. In reality, the biophysiological aspects for neurovascular coupling are far more complicated than the spatiotemporal modulation model in Eq. (1), which deserves a long-term exploration.

Upon a predefined χ source map, we implement 3D BOLD fMRI numerical simulation based on MRI principles. In the results, we are concerned with pattern comparison between the source distribution and the output images (magnitude and phase). Our simulations show that neither the magnitude image nor the phase image is a faithful representation of the χ source distribution. In the context of volumetric medical imaging, the MRI output image is not a tomographic reproduction (quantitative spatial mapping) of the χ source. Since the sourceimage mismatch is due to a cascade of MRI transformations that cause distortions during data acquisition, this inspires us to reconstruct the χ source by solving an inverse MRI problem [21, 22].

poral modulation model in Eq. (1). In particular, our simulation shows the correlation saturation in extreme cases of little noise or noiseless settings; this phenomenon may be explained by the scale invariance of correlation coefficient. On the other extreme case, a severe noise may destroy the task-correlated activity blob; this explains the pursuit on high-SNR

Our 4D BOLD fMRI simulations show that the predefined BOLD NAB in **Figure 13(a)** could be largely reproduced by a task-correlation magnitude fmap (*A*tcorr in **Figure 15(a)**) as extracted from a 4D BOLD fMRI magnitude dataset, thus justifying the BOLD fMRI experiment for brain functional mapping. In comparison, the phase fmap (in *P*tcorr) is spatially dissimilar to the predefined BOLD NAB due to the conspicuous dipole effect in the phase images [6, 11]. Nevertheless, our 4D BOLD fMRI simulations also show that the imaging noise has a strong effect on the fmap extracted from the magnitude or phase image dataset due to the simplified

The data acquisition of BOLD fMRI is not analytically tractable due to the involvement of diversified parameters. The BOLD fMRI simulations provide a means to observe the effect of MRI transformations on the MR data acquisition; spatial distortions between the underlying magnetic source and MR images; and reproducibility of functional activity extraction from a

Since MRI is designed to measure a magnetic field distribution, the BOLD fMRI only measures the χ-expressed BOLD response during a functional activity, a phenotypic numeric expression of a biophysiological brain functional activity in terms of tissue magnetism. It is believed that a functional activity causes IV blood magnetism change in terms of oxygenated and deoxy‐ genated blood magnetic susceptibility change. Therefore, our simulation begins with a configuration of magnetic source by a vasculature-laden VOI with a 3D χ source distribution. Through a spatiotemporal modulation by a predefined local neuroactive blob (numerically NAB(**r**)) and a task paradigm (numerically task(*t*)), we define a dynamic χ source to represent a χ-expressed BOLD activity (in Eq. (1)). It is pointed out in Eq. (1) that the BOLD χ response may incorporate the factors of cerebral metabolic rate of oxygen consumption (CMRO2), cerebral blood volume (CBV), and cerebral blood flow (CBF) through the parameters of *Y*(**t**) and *V*(**r**,*t*), thereby enabling the numeric simulations of MRI-detected BOLD activity. In reality, the biophysiological aspects for neurovascular coupling are far more complicated than the

spatiotemporal modulation model in Eq. (1), which deserves a long-term exploration.

Upon a predefined χ source map, we implement 3D BOLD fMRI numerical simulation based on MRI principles. In the results, we are concerned with pattern comparison between the source distribution and the output images (magnitude and phase). Our simulations show that neither the magnitude image nor the phase image is a faithful representation of the χ source distribution. In the context of volumetric medical imaging, the MRI output image is not a tomographic reproduction (quantitative spatial mapping) of the χ source. Since the source-

spatiotemporal modulation model for numerical BOLD χ expressions.

image acquisition.

22 Numerical Simulation - From Brain Imaging to Turbulent Flows

**4. Discussion**

4D BOLD fMRI dataset.

In NMR principle, a voxel signal is formed from numerous hydrogen proton precessions in a magnetic field. The signal formation involves a huge space scale span from a microscopic atomic scale to macroscopic millimeter scale. For numerical simulations, we implement the mesoscopic micrometer scale through gridel sampling [14, 15]. A gridel is a tiny grid element (at micronmeter scale) smaller than vessel size with which we may digitally depict a vessel geometry. On the other hand, a gridel consists of numerous protons at microscopic atomic scale. The collective proton spins in a gridel are denoted by a spin packet [14, 15]. We define a cortex VOI in a large finely-gridded matrix and partition the VOI into a coarse multivoxel matrix, with each voxel containing an adequate number of gridels for subvoxel structure representation. The VOI partition and gridel rebinning for multivoxel image formation is a topic of multiresolution BOLD signal analysis [15, 23].

In the past decades, the BOLD fMRI mechanism was numerically simulated with signal voxel signals, [7–9], offering an understanding of BOLD fMRI signal formation with respect to a diverse set of parameter settings. However, the single-voxel signal simulation cannot reveal the spatial context for source-image mapping study. Therefore, we were motivated to use multivoxel image simulations for revealing the spatial mismatches between the source and the image [9–11]. Based on 3D BOLD fMRI simulations, it is a straightforward process to imple‐ ment 4D BOLD fMRI simulations. Our 4D BOLD fMRI simulations for a task-evoked brain functional activity, based on a simple spatiotemporal modulation model in Eq. (1), show that the fmap extraction from a 4D BOLD fMRI dataset is sensitive to the additive Gaussian noise. The noise dependence of the task-correlation-based fmap extraction is attributed to the scale invariance of the correlation coefficient.

One factor for the source-image mismatch is the dipole effect that is introduced during the tissue magnetization in a main field *B*0, which is unavoidable for MRI data acquisition. The dipole effect is introduced to the χ-induced fieldmap, which is propagated to the MR magni‐ tude and phase images (signals) via different data transformations. The dipole effect on the χ-induced fieldmap manifests bipolar-valued quadruple lobes around vessels. Upon MRI data acquisition, the magnitude image is a nonnegative nonlinear spatial mapping of the fieldmap and the phase image is an *arctan* nonlinear spatial mapping. It is interesting to show in our numeric simulations that the nonnegative magnitude image resembles the predefined χ source distribution, except for the negative inversions at negative χ regions (not reported herein), and we have found that the (unwrapped) bipolar-valued phase image conforms very well with the bipolar-valued fieldmap [6, 11]. The phase image bears a conspicuous dipole effect that makes a striking difference between the phase image and the predefined χ source.

BOLD fMRI simulation is a time-consuming computation job. In a computer cluster (a Kernel Linux system with 16 CPUs and 252 GB memory), the single-voxel signal simulation requires about 1 hour, and 3D multivoxel simulation requires more than 10 hours, and 4D BOLD fMRI simulation requires a few days, depending on the sizes of gridel, voxel, and VOI. The compu‐ tation burden may be greatly reduced by a Bloch technique [24], which implements intravoxel dephasing signal calculation by a linear approximation. The fast Bloch simulation method is good for MR phase image simulation, but not good for MR magnitude simulation due to an accentuated edge effect. Moreover, the IV and EV signal separation and the diffusion simula‐ tions are not implementable by the Bloch method.
