## 1. Introduction

Magnetic resonance imaging (MRI) is a versatile non-invasive imaging technology that has been widely accepted for brain imaging (probing a magnetic state of brain interior). When applied to brain functional imaging, MRI produces a timeseries of images that are construed as an image representation of a brain functional activity. It is believed that any brain activity incurs a cerebral blood oxygenation level dependent (BOLD) magnetic state change that can be detected by MRI [1–4]. Brain functional imaging based on MRI and the endogenous BOLD contrast is termed BOLD fMRI.

In principle, the MRI output is a complex-valued image consisting of a pair of magnitude and phase [5]. Nevertheless, only the MR magnitude image has been exploited for brain imaging (structural or functional). Recent research shows that neither the MR magnitude nor the phase could faithfully represent the brain magnetic state. This is due to a cascade of MRI transformations (including linear dipole-convolved magnetization and nonlinear complex modulo/argument operations [6]). Consequently, conventional BOLD fMRI that is based on MR magnitude imaging may deviate from the underlying brain magnetic source change due to nonlinear data transformations associated with MR magnitude image formation. Since there is a lack of analytic formulation for describing the imaging aspects of BOLD fMRI, we conduct numeric simulations to understand the BOLD fMRI model.

In the past decades, there have been reports on single-voxel BOLD signal simulations [7–9] and multivoxel 3D BOLD imaging simulations [6, 10, 11]. In this chapter, we first provide a tutorial on the numeric simulations of single-voxel signals and multivoxel images and move forward to address implementing 4D BOLD fMRI simulations.

## 2. Models and methods

An overview of a brain BOLD fMRI model is diagramed in **Figure 1**, which consists of a cascade of three modules (stages). Specifically, the “**Source Magnetism**” module provides the phenotypical χ expression of a brain functional biophysiological state, which serves as the source of the “**MRI technology**” module that produces a complex-valued MR image. Upon data acquisition of a 4D fMRI, a postprocessing stage of “**Statistic image analysis**” is performed to extract the brain functional map (fmap). A complete BOLD fMRI simulation implements the three cascaded stages in **Figure 1** by numerical representations and computations.

### 2.1. Definition of 3D vasculature and magnetic susceptibility source (χ)

The initial step of BOLD fMRI simulation is to configure a χ-expressed BOLD activity, thereby providing a BOLD χ source for fMRI. We define a brain cortex volume of interest (VOI) with a tissue background and fill it with a cortex vasculature, thus simulating a brain cortex region. Let χ_{0}(**r**) denote a static 3D χ distribution of parenchymal tissue in VOI and Δχ(**r**,*t*) the vascular blood χ change associated with a BOLD activity, with **r **= (*x,y,z*) denoting the spatial coordinates in VOI, then the dynamic 4D χ source is given by

(1) |

where* Hct* denotes the blood hemocrit (*Hct* = 0.4 for normal blood), χ_{do} the magnetic susceptibility difference between deoxygenated and oxygenated blood tissues (χ_{do} = 0.27 × 4π ppm (in SI unit)), *Y*(*t*) the blood oxygenation level (*Y* ∈ [0,1]), NAB(**r**) the local neuroactive blob distribution, and *V*(**r**,*t*) the vasculature geometry in VOI. The explicit *t* variable indicates a possible change during a BOLD activity, such as cerebral blood volume change in *V*(**r**,*t*) and oxygenation level change in *Y*(*t*). For the sake of simulating fMRI signals, a pure BOLD activity is expressed by a dynamic blood magnetic susceptibility change, Δχ(**r**,*t*), which serves as the magnetic source for BOLD MRI simulation. In practice, the BOLD activity provides an additive term, Δχ(**r**,*t*) (a perturbation term), on a background distribution χ_{0}(**r**) in Eq. (1).

A local functional activity is defined by a 3D spatial distribution of NAB(**r**) (a neuroactive blob centered at **r** in VOI). For the sake of numerical representation, we assume a NAB by a Gaussian-shaped blob (with soft boundary) or a ball-shaped blob (with hard boundary). A NAB defines a local activity distribution in VOI, which presents with an ON state (active state) and vanishes with an OFF state (resting state) by a temporal modulation of a designed task paradigm. We may define an excitatory activity by a positive distribution (NAB(**r**) > 0) or an inhibitory activity by a negative distribution (NAB(**r**) < 0), in relation to the static background distribution. For the numerical simulation of a BOLD activity, we define a BOLD χ response by a spatiotemporal modulation model in Eq. (1). A brain active state gives rise to Δχ(**r**,*t*) ≠ 0 in NAB and at a task “ON” epoch, and a brain resting state is numerically characterized by Δχ(**r**,*t*) = 0 over the VOI in Eq. (1).

It is mentioned that the BOLD χ expression in a brain activity is simply simulated by a spatial modulation model in Eq. (1), where a neuronal activity is defined by a local blob that shapes a local blood Δχ map by a spatial multiplication. We also simplify the BOLD χ source simulation by ignoring the hemodynamic response function (hrf), which otherwise could be accounted for by convoluting Δχ(**r**,*t*) with a kernel of hrf (usually adopting a canonical hrf that is characteristic of a high upshoot followed by a small undershoot).

A BOLD χ change happens inside the vascular blood stream. We need to configure the vasculature geometry *V*(**r**,*t*) by filling the VOI with cluttered vessels with a blood volume fraction (*bfrac*), as expressed by

where the *t* variable is reserved to incorporate the change in cerebral blood volume as a result of vasodilation/vasoconstriction in a BOLD activity. A static vasculature is included as a binary volume *V*(**r**) that remains stationary during a BOLD activity. The random vascular geometry is generated under a control of *bfrac *= [0.02, 0.04] for cortex vasculature simulation [1, 8, 11–13].

In order to maintain a control of constant *bfrac* for cortex vasculature over different regions or across multiresolution subregions, we may fill a VOI with random beads instead of cluttered vessels. In **Figure 2** (**a1**,**b1**) are illustrated two brain local vasculature geometries with cluttered cylinders and random beads, under local stimuli by an excitatory blob (in red) and an inhibitory blob (in green). The NAB-modulated BOLD χ response distributions (in an active ON state) are shown in **Figure 2** **(a2, b2)** with a *y*_{0}-slice in which the inactive regions (far from NAB) have little or no BOLD responses (Δχ ≈ 0).

In order to numerically depict the vasculature geometry, we need to define the VOI with a large finely gridded 3D matrix with a tiny grid element (gridel) at a scale of micronmeter [14]. For example, a matrix of 2048 × 2048 × 2048 gridels, where a gridel = 2 × 2 × 2 μm^{3}, is used to represent a small VOI of 4.1 × 4.1 × 4.1 mm^{3}. The large matrix resulting from VOI gridel sampling offers a quasi-continuous representation of a continuous distribution over VOI. A gridel represents a spin packet (or isochromat) that contains numerous identical proton spins, serving as a mesoscopic representation (at micronmeter scale) between microscopic structure (at atomic and molecular angstrom scale) and macroscopic structure (at millimeter scale of MRI voxels) [15, 16].

### 2.2. Calculation of χ-induced fieldmap

Upon determining the brain χ source configuration, we calculate the χ-induced magnetic field map (fieldmap for short) by a 3D spatial convolution with a dipole kernel. This is to simulate the brain tissue magnetization process in a main field B_{0} that produces an inhomogeneous fieldmap. Let b(*x,y,z*) represent the z-component of χ-induced 3D vector field; it is given by [10, 11]

where * denotes spatial convolution, and *h*_{dipole} a 3D dipole field [17]. In a Fourier domain, the 3D dipole convolution can be efficiently implemented by multiplicative spatial filtering, as expressed by [18]

(4) |

where (*k _{x}*,

*k*,

_{y}*k*) denotes the coordinates in the Fourier domain. The fieldmaps induced by the Δχ distribution in a main field

_{z}*B*

_{0}are illustrated in

**Figure 3**(displayed with a

*y*

_{0}-slice), which shows a conspicuous dipole effect in a manifestation of bipolar-valued quadruple lobes around vessels (with an orientation dependence [19]).

In computation implementation, the 3D FFT for fieldmap calculation for a finely-gridded 3D χ distribution matrix (e.g., 2048 × 2048 × 2048 gridels) may encounter an “out-of-memory” problem. We propose to solve this problem by decomposing 3D FFT into 2D FFT and 1D FFT. Specifically, we first conduct 2D FFT on each *z*-slice (or *xy*-plane) and save the data as data files, and then conduct 1D FFT along each of the *z*-axis columns of a 3D volume that is stacked from *z*-slices (processed by 2D FFT and saved in files). In order to reduce the data file management (*fwrite* and *fread*) of the 3D FFT decomposition, we decompose the 3D matrix into a handful of *z*-chunks (z-slabs) that each consists of multiple *z*-slices. The number of *z*-chunks is dependent upon the available computer RAM (random access memory). As illustrated in **Figure 4**, we only need to manage (*fwrite *and *fread*) a number of 32 *z*-chucks (each consists of 64 *z*-slices in a matrix of 2048 × 2048 × 64), instead of 2048 individual *z*-slices.

### 2.3. Multivoxel partition of VOI

An MRI output is a discrete multivoxel image with the voxel size at a macroscopic millimeter scale, which implies that the MRI scanning process partitions a brain VOI into a small array of macroscopic voxels. We simulate a multivoxel MR image by rebinning mesoscopic gridels (at micronmeter scale) into macroscopic voxels (at millimeter scale). For example, we can reduce a large matrix of 2048 × 2048 × 2048 gridels to a small image matrix of 64 × 64 × 64 voxels with a voxel of 32 × 32 × 32 gridels. The multivoxel partition of VOI is in fact a spatial sampling by voxels, called voxelization. We denote the gridel-sampled representation of a continuous distribution over VOI by a spatial variable “(**r**)”, and the voxel-sampled discrete representation by an index variable “[**r**]”. Let Ω denote a voxel space, and |Ω| the voxel size (in terms of a number of gridels in Ω). The VOI voxelization is represented by

The VOI voxelization (voxel sampling) is necessary for MRI to produce a multivoxel image, due to the band limit of coil transmission and reception, which is designed as a parameter of voxel size in MRI protocol. The voxel size also represents a parameter of spatial resolution. In the MRI output, a high-resolution (corresponding to small voxel size) produces a large image matrix, and vice versa.

### 2.4. Calculation of intravoxel dephasing signals (Monte Carlo method)

An MRI voxel signal (or a NMR signal) is formed by an intravoxel spin precession dephasing in a χ-induced fieldmap. A quadrature detection produces a complex-valued voxel signal, denoted by C[*x,y,z*] that is formulated by [5]

where γ denotes the gyromagnetic ratio, and the auxiliary variable ‘X’ is reserved to explicitly include the dependence of NMR signal upon a diverse set of factors. We are always concerned with the factors of echo time (*T _{E}*), field strength (

*B*

_{0}), spatial resolution (voxel size |Ω|), and vessel geometry in particular.

A voxel contains a number of gridels that each represents a spin packet. The voxel signal calculation in Eq. (6) counts all the spin packets in the voxel space. For a voxel that contains a large number of gridels, we may select a smaller number of gridels to calculate the voxel signal and reduce the computation burden. The intravoxel dephasing signal calculation made by counting the spin packets is a Monte Carlo simulation, which is expressed by

For example, a voxel of 32 × 32 × 32 gridels consists of 32,768 spin packets, from which we may randomly select 3000 for the intravoxel average computation in Eq. (7) at a small computation cost of 10% (≈3000/32,768). It is noted that *C*[*x,y,*z;*X*] denotes a complex voxel signal at [*x,y,z*] in VOI, and we also use *C*[*x,y,z*;*X*] to represent a multivoxel complex-valued image in the context that [*x,y,z*] addresses all the voxels in VOI.

From a complex signal (image), we can calculate its magnitude and phase components by

It is also noted that we use the magnitude loss and phase accrual to represent the pair of complex signal components and that the magnitude and phase calculations are different nonlinear operations.

### 2.5. Intravascular (IV) and extravascular (EV) signal separation

In an MRI experiment, it is difficult to separate intravascular (IV) signal from extravascular (EV) signal in an NMR signal. In numerical simulation, we can calculate the IV and EV signals separately based on the binary partition of voxel space according to the vessel geometry. Let Ω^{IV} and Ω^{EV} denote the IV and EV subspaces in a voxel space, which are partitioned by the vessel geometry by

(9) |

Then, we calculate the IV signal by only counting the gridels that are within vessel space (Ω^{IV}), and the EV signal by the gridels in Ω^{EV}. That is, the IV and EV signals are given by

(10) |

In **Figure 5** are illustrated the IV/EV partition of a voxel space for IV/EV signal simulations. It is mentioned that the Ω^{IV} only occupies a small fraction of Ω and the BOLD χ change is confined in Ω^{IV}.

Although a BOLD Δχ change is confined in Ω^{IV} in a NAB, the vascular blood magnetization process in B_{0} establishes a long-range magnetic field distribution, not only in Ω^{IV} but also in Ω^{EV}, with a distant decay (∝1/*r*^{3}) and a spatial modulation by NAB (see **Figures 2** and **3**). Obviously, a BOLD activity causes an IV signal and an EV signal simultaneously, which are generated from different field values over the IV and EV spaces, respectively. In **Figure 5(a)** is illustrated the IV/EV signal formations from spin particles in the IV/EV partition spaces.

A voxel NMR signal is formed from its IV and EV signals by a convex combination according to the IV/EV occupancies, as represented by

Consequently, the IV signal contribution is greatly suppressed by a small weight of *bfrac* (<<1), as will be demonstrated later.

### 2.6. Diffusion effect

An NMR signal is formed via the carrier of hydrogen protons in tissue water. Since the water molecules undergo random motions, the water protons are non-stationary. We describe the proton random motion in 3D space by a trajectory **r**(*t*), which is represented by [9, 20]

(12) |

where *δt* denotes the time interval of the random motion of water molecules (*δt* = 2 ms in simulation), *D* the diffusion coefficient (different for diffusions in IV and EV), and Gauss a Gaussian distribution of the random motions (with a standard deviation of *σ _{d}*). It is noted that water proton diffusion in IV space is twice faster than in EV space. In

**Figure 5(b)**are illustrated the diffusion IV and EV signal simulations.

### 2.7. Volumetric BOLD fMRI simulation

Based on individual voxel signal calculation, we implement 3D volumetric BOLD fMRI simulation by calculating the voxel signals at a multivoxel image array. Given a 3D χ source, the 3D BOLD fMRI simulations produce a 3D complex-valued multivoxel image *C*[*x,y,z*; *X*]; here, we are concerned with the spatial pattern comparison between the χ source and the magnitude image. Since the phase image bears a conspicuous dipole effect that dooms the morphological mismatch with the χ source, we do not need to compare the phase image with the χ source. However, the phase image is directly related to the χ-induced fieldmap, and the phase image has been used for the fieldmap reconstruction in an inverse MRI solver [11, 21, 22]. In particular, in a small phase angle regime, the phase conforms the fieldmap with a difference of constant scale. In large phase angle scenarios, the unwrapped phase image resembles the fieldmap very well (albeit with somewhat nonlinear distortions). Therefore, we are also concerned with pattern comparison between MR phase image and the fieldmap. We suggest the spatial pattern comparisons by spatial correlations, which are defined by

It is noted that the spatial pattern correlations are applied to the multivoxel matrices (in notation of [**r**]) of the χ source, the fieldmap, and the images, which are all discretized at a spatial resolution of the same voxel size.

### 2.8. 4D BOLD fMRI simulation

It is straightforward to implement 4D BOLD fMRI simulation based on the repetition of volumetric BOLD fMRI simulations for each snapshot capture over a BOLD activity. First, we need to define a task-evoked 4D BOLD χ change, as illustrated in **Figure 6**. Specifically, we configure a 3D vasculature-laden VOI and provide a 3D χ distribution for a brain VOI state. A local χ change is simulated with a spatiotemporal modulation by a NAB and a task paradigm (in Eq. (1)). For the weak BOLD response detection, the task paradigm is usually designed as a boxcar waveform for repetition measurement of BOLD signals. We may define a positive NAB for an excitatory BOLD response and a negative NAB for an inhibitory response. The static background χ_{0} may be assigned to a water pool (χ_{0} = χ_{water}) or be empty (χ_{0} = 0) with an additive Gaussian noise.

The 4D BOLD fMRI simulation involves a predefined 4D source χ[**r**,*t*] and two output 4D images (*A*[**r**,*t*] for magnitude and *P*[**r**,*t*] for phase, as defined in Eq. (8)). Conventional BOLD fMRI exploits the 4D magnitude dataset *A*[**r**,*t*] for functional analysis. For a task-evoked BOLD fMRI simulation, the functional activity map can be extracted from a 4D dataspace, Λ[**r**,*t*] = {χ[**r**,*t*], *A*[**r**,*t*], *P*[**r**,*t*]} by a temporal correlation (*tcorr)* map that is defined by

(14) |

where std_{t} denotes the standard deviation of the data with respect to the *t* variable, χ^{true} the predefined χ source, and χ^{recon} the reconstructed χ source (by solving the inverse problem of MRI data). It is noted that the correlation coefficient is invariant to signal strength. Therefore, a strong response signal may have the same correlation value as a weak response does as long as the strong and weak responses take on the same timecourse profiles. Consequently, the scale invariance of correlation leads to correlation saturation (*tcorr* = 1 at regions with different response strengths). Nevertheless, the correlation saturation can be ruined by the presence of noise. Herein, by noise we mean any pattern difference between the response signal timecourse and task timecourse. In reality, the BOLD χ responses are subject to various noises (biological noise, physiological noise, detection noise) that spoil the task correlations at weak response regions. Only strong responses are immune to noise spoilage. It is the noise in the voxel response timecourse (extracted from a 4D dataset at a specific voxel) that shapes the *tcorr* map according to the response signal strength.

## 3. Simulation results

### 3.1. Single voxel signal simulations

#### 3.1.1. EV/IV signal separation

By calculating the EV and IV signal portions separately and their convex linear combination in Eq. (11), we present the EV/IV signal behaviors with respect to a span of echo time (*T _{E}*=[0, 60] ms) and for a range of field strength (

*B*

_{0}=[1.5, 3, 4, 7, 9] Tesla). It is seen in

**Figure 7**that the IV signal changes quickly with a long echo time. However, the drastic IV signal changes are greatly suppressed in the voxel signals by the dominant EV signals. In particular, the IV signal may be developed into phase wrapping phenomenon for a long echo time (see

**Figure 7(b2)**). With the dominance of EV signals in a large voxel, a voxel signal remains as a linear phase accrual with echo time (see

**Figure 7(b3)**).

#### 3.1.2. Multiresolution voxel signal behavior

As a voxel size decreases, the voxel space contains less (or none) vessels, and there is less voxel average effect. In **Figure 8**, the four-level voxel subdivision and multiresolution voxel signal behaviors are demonstrated. At level =1, the parent voxel contains a clutter of vessels where the complex voxel signal appears as a short line-segment trajectory (with respect to *T _{E}*). As the voxel is decomposed into an 8 × 8 × 8 array at level = 4, the subvoxel only contains a single vessel, and the voxel signal becomes turbulent due to the high field values for rapid Larmor precession [14, 23].

#### 3.1.3. Diffusion effects on magnitude and phase signals

The numerical simulations on the diffusion effect on MR magnitude and phase are presented in **Figure 9** for a span of *T _{E}* = [0, 60] ms with different field strengths (in terms of Δχ

*B*

_{0 }= [0.1, 3] ppmT). The results show that the diffusion has more effect on low field magnitude than on high field magnitude [20]. Nevertheless, the diffusion has little effect on MR phase signals.

### 3.2. Volumetric BOLD fMRI simulations

#### 3.2.1. Cortex VOI configuration and voxelization

We define a cortex VOI in a large matrix and fill it with random beads (radius = 3 μm, *bfrac *= 0.03), and simulate local BOLD response by a Gaussian-shaped NAB, which modulates the local χ distribution by a spatial multiplication. The VOI is partitioned into a coarse matrix by grouping the gridels into voxels. As a result of voxelization, we represent a distribution over VOI by a multivoxel image matrix. The voxelization with a large voxel size produces a small image matrix, and vice versa. **Figure 10** shows a VOI that is represented by a large matrix in gridel sampling **(a)** with a zoomed region for substructure visualization **(b)**. The VOI voxelization by a voxel of 32 × 32 × 32 gridels produces a matrix of 64 × 64 × 64 voxels **(c)** and produces a matrix of 32 × 32 × 32 voxels **(d)** by a voxel of 64 × 64 × 64 gridels. It is seen that a larger voxel size is comprised of more spatial smoothing.

#### 3.2.2. Multivoxel image calculation

Given a 3D χ distribution in **Figure 11(a)**, we calculated the χ-induced fieldmap by Eq. (2) and presented the results **(b)**. In the absence of diffusion, we calculated the complex-valued T2* images (**c**, **d**). In the presence of diffusion (Eq. (12)), we recalculated the complex-valued T2* images (**e**, **f**). The diffusion simulation on multivoxel fMRI shows that the diffusion effect is insignificant on image formation.

#### 3.2.3. Morphological distortions associated with 3D BOLD fMRI

We performed volumetric BOLD fMRI simulations for a span of echo times (*T _{E}*

_{ }= [0, 30] ms) with different parameter settings with respect to voxel size, field strength, and with and without diffusion. With the datasets of numerical BOLD fMRI simulations, we compared the magnitude images with the predefined χ source and the phase images with the fieldmaps. The results are presented in

**Figure 12**. Note that the pattern correlations are plotted in a small display window ([0.9, 1]) out of the full range of

*corr*∈ [–1, 1] for scrutiny.

### 3.3. 4D BOLD fMRI simulations

The 4D BOLD fMRI simulations are presented in **Figures 13** through **15**. Specifically, in **Figure 13** are shown **(a)** the VOI configuration with two local neuroactive blobs (NAB), **(b)** the NAB-modulated BOLD χ distribution at an ON state (or active state), and **(c)** the NAB-absent χ distribution at an OFF state (or resting state), displayed with a *y*_{0}-slcie. We designed a task paradigm by a pattern of 5 ON states and 5 OFF states, simulating the brain active state measurement by 5 repetitions and the brain resting state measurement by another 5 repetitions. (In practice, a multiple repetition of the “ON/OFF” pattern is used to design the task paradigm). The bead-represented vasculature structure in a voxel in VOI is shown in zoom **(d)** with a 3D display. It is noted that the VOI is represented in a matrix of 2048 × 2048 × 2048 gridels **(a)**, the voxelization on VOI is represented by a multivoxel matrix of 64 × 64 × 64 voxels **(b)** and **(c)** with a voxel = 32 × 32 × 32 gridels **(d)**, and that the cortex vasculature in a VOI is simulated by a uniform random distribution (background) that is independent of the BOLD NAB and task paradigm. The 4D χ(**r**,*t*) representation for a local BOLD activity is related to the NAB and the task through a spatiotemporal model in Eq. (1).

Upon the numerical representation of 4D χ(**r**,*t*), we performed 4D BOLD fMRI simulations by repeating the 3D BOLD fMRI simulation for each snapshot time point (there are 10 timepoints for the task pattern of 5 ONs and 5 OFFs), with the settings (*T _{E} *= 30 ms,

*B*

_{0}= 3 T, VOI matrix = 64 × 64 × 64 voxels, 1 voxel = 32 × 32 × 32 gridels, 1 gridel = 2 × 2 × 2 μm

^{3}).

**Figure 14**shows (

**a1**,

**a2**) the magnitude images, (

**b1**,

**b2**) the phase images captured at an ON and OFF state, and (

**a3**,

**b3**) the timecourses of magnitude signal changes, and phase signal changes at two voxels: one voxel inside an active blob (marked by “x” in (

**a1**)) and another outside the active at blob (marked by “o”). It is noted the ripples in the signal timecourses in (

**a3**,

**b3**) are attributed to the additive Gaussian random noise in the data acquisition simulations.

By arranging the timeseries of images according to the task timecourse, we can play a movie for a BOLD activity. In reality, the BOLD response is too weak and noisy to be perceived between an ON and OFF state. For the sake of BOLD response pattern representation, we need to extract the BOLD activity blobs from the timeseries of images by statistical parameter mapping method, which consists of an essential task correlation map as defined in Eq. (14).

Upon the completion of 4D BOLD magnitude and phase image datasets (*A*[**r**,*t*], *P*[**r**,*t*]), we calculated the task- correlated fmap using Eq. (14). In the results, we obtained 3D *A*_{tcorr} for BOLD magnitude fmap from the 4D magnitude image dataspace, and a 3D *P*_{tcorr} for BOLD phase fmap from the 4D phase image dataspace. By repeating the 4D BOLD fMRI simulations with different noise levels, we show that *A*_{tcorr} or *P*_{tcorr} is sensitive to the additive Gaussian noise. In **Figure 15** are showed the *A*_{tcorr} and *P*_{tcorr} (displayed with a *y*_{0}-slice out of the 64 × 64 × 64 matrix volume) for the Gaussian noise at different noise levels = {0.001, 0.01, 0.05, 0.1}. It is seen that either the magnitude or phase fmap suffers from correlation saturation in little noise (noise level < 0.01) or tends to be buried in severe noises (noise level > 0.05) for our spatiotemporal 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 image acquisition.

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 spatiotemporal modulation model for numerical BOLD χ expressions.

## 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 deoxygenated 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 (CMRO_{2}), 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-image 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].

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 implement 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 magnitude 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 computation 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 simulations are not implementable by the Bloch method.

## 5. Conclusion

We conclude our numerical BOLD fMRI simulations by the following findings (albeit qualitative):

Both the MR magnitude and the phase images are spatially different from the predefined magnetic susceptibility distribution. This image-source distortion is due to the inevitable data transformations associated with MRI data.

By numerical simulation, we can separate the intravascular (IV) signal from the extravascular (EV) signal in a voxel signal. The IV signal is much stronger than the EV signal as a result of a BOLD χ change. However, the drastic IV signal evolution is usually greatly suppressed in a voxel signal by a small proportion of blood volume fraction (

*bfrac*≈ [0.02, 0.04]).As voxel size decreases, the voxel signals evolve more drastically and turbulently inside and around the large vessels.

The proton diffusion effect due to nonstationary water molecules in brain tissues incurs more MR magnitude signal decays in a low field than in a high field. In comparison, the proton diffusion has little effect on MR phase signals.

The numerical simulation on 4D BOLD fMRI for task-evoked functional mapping shows that the functional activity extraction by a task correlation technique is sensitive to data noise.

Overall, the numerical simulations on BOLD fMRI allow us to look into the insights of a single-voxel signal, a multivoxel image, and a video of brain functional BOLD activity with respect to various parameter settings. The finding in source-image mismatch inspires us to seek for the underlying magnetic source of BOLD fMRI for more accurate brain functional mapping. The finding in the noise sensitiveness of task-correlated fmap raises a caveat to the correlation-based functional mapping.