**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 B0 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]

$$\begin{aligned} b(\mathbf{x}, \mathbf{y}, z) &= B\_0 \chi(\mathbf{x}, \mathbf{y}, z) \ast h\_{\text{dipole}}(\mathbf{x}, \mathbf{y}, z) \\ \text{with } \quad h\_{\text{dipole}}(\mathbf{x}, \mathbf{y}, z) &= \frac{1}{4\pi} \frac{3z^2 - \left(\mathbf{x}^2 + \mathbf{y}^2 + z^2\right)^{3/2}}{\left(\mathbf{x}^2 + \mathbf{y}^2 + z^2\right)^{5/2}} \end{aligned} \tag{3}$$

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]

$$\begin{aligned} FT\{b(\mathbf{x},\mathbf{y},\mathbf{z}) &= H\_{\text{dipole}}(k\_x, k\_y, k\_z) \cdot FT\{\boldsymbol{\chi}(\mathbf{x},\mathbf{y},\mathbf{z})\} \cdot B\_0 \\ \text{with } H\_{\text{dipole}}(k\_x, k\_y, k\_z) &= FT\{h\_{\text{dipole}}(\mathbf{x},\mathbf{y},\mathbf{z})\} = \frac{1}{3} - \frac{k\_z^2}{k\_x^2 + k\_y^2 + k\_z^2} \end{aligned} \tag{4}$$

where (*kx*,*ky*,*kz*) denotes the coordinates in the Fourier domain. The fieldmaps induced by the Δχ distribution in a main field *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 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

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].

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

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 B0 that produces an inhomogeneous fieldmap. Let b(*x,y,z*) represent the z-component of χ-induced 3D vector field; it is given by

2 2 2 2 3/2

2

*z*

*xyz*

*xyz*


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

dipole 0

c

º = - + +

dipole dipole 222

where (*kx*,*ky*,*kz*) denotes the coordinates in the Fourier domain. The fieldmaps induced by the Δχ distribution in a main field *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

<sup>1</sup> with ( , , ) { ( , , )} 3

= ××

*<sup>k</sup> H k k k FT h x y z kkk*

0 dipole

( , , ) ( , , ) ( , , )

*bxyz B xyz h xyz*

{ ( , , } ( , , ) { ( , , )}

*FT b x y z H k k k FT x y z B*

*xyz*

vessels (with an orientation dependence [19]).

*xyz*

= \*

c

dipole 2 2 2 5/2

13 ( ) with ( , , ) 4( )

p

*z xyz h xyz*

, is used to

(3)

(4)

. The large matrix resulting from VOI gridel

For example, a matrix of 2048 × 2048 × 2048 gridels, where a gridel = 2 × 2 × 2 μm3

little or no BOLD responses (Δχ ≈ 0).

MRI voxels) [15, 16].

[10, 11]

expressed by [18]

represent a small VOI of 4.1 × 4.1 × 4.1 mm3

8 Numerical Simulation - From Brain Imaging to Turbulent Flows

**2.2. Calculation of χ-induced fieldmap**

**Figure 3.** The fieldmaps calculated from the Δχ distributions in **Figure 2(a2, b2)**. It is noted that the Δχ-induced field‐ map takes on a continuous inhomogenous bipolar-valued distribution over VOI, bearing a conspicuous dipole effect around large vessels (beads).

**Figure 4.** 3D FFT implementation by 2D FFT and 1D FFT. The 3D FFT of a large 3D matrix (e.g., 2048 × 2048 × 2048) is achieved by first performing 2D FFT on each *z*-slice (*xy*-plane) and then 1D FFT along *z* columns. A large 3D matrix is decomposed into a number of small *z*-chunks to reduce the data file management (*fwrite* and *fread*).

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 man‐ agement (*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

$$V[\mathbf{x}, \mathbf{y}, \mathbf{z}] = \frac{1}{|\Omega|} \sum\_{(\mathbf{x'}, \mathbf{y'}, \mathbf{z'}) \in \Omega(\mathbf{x}, \mathbf{y}, \mathbf{z})} V(\mathbf{x'}, \mathbf{y'}, \mathbf{z'}) \tag{5}$$

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]

$$\begin{aligned} \lfloor C[\mathbf{x}, \mathbf{y}, z; X] = \frac{1}{|\Omega|} \sum\_{(\mathbf{x}^\*, \mathbf{y}^\*, z^\*) \in \Omega(x\_\*, z)} e^{\langle \mathbf{y}^\*(\mathbf{x}^\*, \mathbf{y}^\*, z^\*)T\_\mathbf{x} \rangle} \\ \text{with} \quad X = \{ \text{"\$T\_E\$", "B\_0\$", "|\Omega|} \text{"}, \text{"\$\infty\text{-}\$} \text{else} \text{else} ", \cdots \text{]} \end{aligned} \tag{6}$$

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 (*TE*), 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

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

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

> ( , , ) (,,) <sup>1</sup> [, ,] (, ,) | | *x y z xyz Vxyz Vxyz* ¢¢¢ ÎW

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

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,

( ', '. ') ( , , )

*x y z xyz*

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 (*TE*), field strength (*B*0), spatial resolution (voxel size |Ω|), and

å

ÎW

0

<sup>=</sup> <sup>W</sup> = W

*CxyzX e*

*E*

with {' ',' ',' | | ',' ', }

*X T B vesselsize*

( ', ', ')

*<sup>E</sup> ibx y z T*

L

(6)

g

<sup>=</sup> ¢¢¢ <sup>W</sup> å (5)

64 *z*-slices in a matrix of 2048 × 2048 × 64), instead of 2048 individual *z*-slices.

number of gridels in Ω). The VOI voxelization is represented by

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

<sup>1</sup> [, ,; ] | |

**2.3. Multivoxel partition of VOI**

10 Numerical Simulation - From Brain Imaging to Turbulent Flows

matrix, and vice versa.

vessel geometry in particular.

denoted by C[*x,y,z*] that is formulated by [5]

$$\mathbb{E}[\mathbf{x}, \mathbf{y}, z; X] = \frac{1}{N} \sum\_{s=1}^{N \cdot |\Omega|} e^{\langle \gamma b(x\_s, y\_s, z\_s) \rangle\_{\mathbb{E}}} \text{ for } (x\_s, y\_s, z\_s) \in \Omega(\mathbf{x}, \mathbf{y}, z) \tag{7}$$

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

$$\begin{cases} \mathcal{A}[\mathbf{x}, \mathbf{y}, z; X] = \mathbb{I} - \left| C[\mathbf{x}, \mathbf{y}, z; X] \right| & \text{(magnitude loss)}\\ P[\mathbf{x}, \mathbf{y}, z; X] = \mathcal{L}C[\mathbf{x}, \mathbf{y}, z; X] & \text{(phase accrual)} \end{cases} \tag{8}$$

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

$$\begin{aligned} \Omega^{\mathcal{W}}(\mathbf{x}, \mathbf{y}, z) &= \begin{cases} \mathbf{l} & (\mathbf{x}, \mathbf{y}, z) \in \mathsf{vessel} \\ \mathbf{0} & \text{otherwise} \end{cases} \\ \Omega^{E^{\mathcal{V}}}(\mathbf{x}, \mathbf{y}, z) &= \begin{cases} \mathbf{l} & (\mathbf{x}, \mathbf{y}, z) \notin \mathsf{vessel} \\ \mathbf{0} & \text{otherwise} \end{cases} = \mathbf{l} - \Omega^{\mathcal{W}}(\mathbf{x}, \mathbf{y}, z) \\ \text{with} \quad \Omega = \Omega^{\mathcal{W}} \cup \Omega^{E^{\mathcal{V}}} \text{ and } |\Omega| = |\Omega^{\mathcal{W}}| + |\Omega^{E^{\mathcal{V}}}| \end{aligned} \tag{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

$$\begin{split} \left| \boldsymbol{C}^{\mathcal{W}} [\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}; \boldsymbol{X}] = \frac{1}{|\boldsymbol{\Omega}^{\mathcal{W}}|} \sum\_{(\boldsymbol{x}^{\prime}, \boldsymbol{y}^{\prime}, \boldsymbol{z}^{\prime}) \in \boldsymbol{\Omega}^{\mathcal{V}}(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z})} \boldsymbol{e}^{\prime \boldsymbol{y} \Phi(\boldsymbol{x}^{\prime}, \boldsymbol{y}^{\prime}, \boldsymbol{z}^{\prime}) \mathcal{T}\_{\mathcal{E}}} \\ \boldsymbol{C}^{\mathcal{E}^{\mathcal{V}}}[\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}; \boldsymbol{X}] = \frac{1}{|\boldsymbol{\Omega}^{\mathcal{E}^{\mathcal{V}}}|} \sum\_{(\boldsymbol{x}^{\prime}, \boldsymbol{y}^{\prime}, \boldsymbol{z}) \in \boldsymbol{\Omega}^{\mathcal{V}}(\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z})} \boldsymbol{e}^{\prime \boldsymbol{y} \Phi(\boldsymbol{x}^{\prime}, \boldsymbol{y}^{\prime}, \boldsymbol{z}) \mathcal{T}\_{\mathcal{E}}} \\ \text{with} \quad \boldsymbol{\Omega} = \boldsymbol{\Omega}^{\mathcal{W}} \cup \boldsymbol{\Omega}^{\mathcal{E}^{\mathcal{V}}} \end{split} \tag{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 B0 establishes a long-range magnetic field distribution, not only in ΩIV but also in ΩEV, with a distant decay (∝1/*r*<sup>3</sup> ) 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

$$\begin{aligned} \text{C[x, y, z; X]} &= b \text{frac} \cdot C^{\prime\prime} \left[ \mathbf{x}, y, z; X \right] + (\mathbf{l} - b \text{frac}) \cdot C^{E\prime} \left[ \mathbf{x}, y, z; X \right] \\ \text{with} \quad b \text{frac} &= \frac{|\Omega^{\prime\prime}|}{|\Omega|} \end{aligned} \tag{11}$$

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

**Figure 5.** Illustration of extravascular (EV) and intravascular (IV) space partition in a voxel for intravoxel spin dephas‐ ing signal simulations (a) in absence of spin diffusions (static spins) and (b) in the presence of spin diffusions. A voxel space is partitioned by its intravoxel binary vasculature into IV (vessel=1) and EV (vessel=0) subspaces.

### **2.6. Diffusion effect**

( ', ', ')

g

g

*E*

*E*

) and a spatial modulation by NAB (see **Figures 2** and **3**).

å (10)

(11)

( ', ', ')

( ', '. ') ( , , )

ÎW

*x y z xyz*

å

*IV ibx y z T IV*

*EV ibx y z T EV*

*IV*

<sup>1</sup> [, ,; ] | |

<sup>1</sup> [, ,; ] | |

W=W ÈW

*IV EV*

<sup>=</sup> <sup>W</sup>

*C xyzX e*

<sup>=</sup> <sup>W</sup>

*C xyzX e*

with

12 Numerical Simulation - From Brain Imaging to Turbulent Flows

ΩEV, with a distant decay (∝1/*r*<sup>3</sup>

as will be demonstrated later.

to the IV/EV occupancies, as represented by


*bfrac*


W = W

*IV*

in ΩIV.

( ', '. ') ( , , )

ÎW

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

Although a BOLD Δχ change is confined in ΩIV in a NAB, the vascular blood magnetization process in B0 establishes a long-range magnetic field distribution, not only in ΩIV but also in

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

Consequently, the IV signal contribution is greatly suppressed by a small weight of *bfrac* (<<1),

**Figure 5.** Illustration of extravascular (EV) and intravascular (IV) space partition in a voxel for intravoxel spin dephas‐ ing signal simulations (a) in absence of spin diffusions (static spins) and (b) in the presence of spin diffusions. A voxel

space is partitioned by its intravoxel binary vasculature into IV (vessel=1) and EV (vessel=0) subspaces.

*IV EV*

[ , , ; ] [ , , ; ] (1 ) [ , , ; ]

*C x y z X bfrac C x y z X bfrac C x y z X*

= × +- ×

*x y z xyz*

*EV*

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]

$$\begin{aligned} \mathbf{r}(t) &= \mathbf{x}(t)\hat{\mathbf{x}} + \mathbf{y}(t)\hat{\mathbf{y}} + \mathbf{z}(t)\hat{\mathbf{z}} \\ \mathbf{x}(t + \delta t) &= \text{Gauss}(\mathbf{x}(t), \sigma\_{d}) \\ \mathbf{y}(t + \delta t) &= \text{Gauss}(\mathbf{y}(t), \sigma\_{d}) \\ \mathbf{z}(t + \delta t) &= \text{Gauss}(\mathbf{z}(t), \sigma\_{d}) \\ \text{with} \quad \sigma\_{d} &= \sqrt{2 \cdot D \cdot \delta t} \\ D &= \begin{cases} 1.5 \times 10^{-3} cm^{2} \text{ / s} & (\mathbf{x}, \mathbf{y}, \mathbf{z}) \in \text{vessel} \\ 0.75 \times 10^{-5} cm^{2} \text{ / s} & (\mathbf{x}, \mathbf{y}, \mathbf{z}) \notin \text{vessel} \end{cases} \end{aligned} \tag{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

$$\begin{aligned} \text{corr}(A, \mathbb{X}) &= \frac{\text{cov}(A[\mathbf{r}], \mathbb{X}[\mathbf{r}])}{\text{std}(A[\mathbf{r}]) \cdot \text{std}(\mathbb{X}[\mathbf{r}])} \\ \text{corr}(P, b) &= \frac{\text{cov}(P[\mathbf{r}], b[\mathbf{r}])}{\text{std}(P[\mathbf{r}]) \cdot \text{std}(b[\mathbf{r}])} \end{aligned} \tag{13}$$

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.

**Figure 6.** Illustration of 4D BOLD χ response simulations. A VOI is filled local Δχ change with positive and negative Δχ responses superimposed on a static background in the presence of noise. A BOLD event is represented by a times‐ eries of the 3D Δχ snapshot distributions in Eq. (1) through a spatiotemporal modulation by NAB(**r**) and task(*t*).

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

$$\Lambda\_{\text{tower}}[\mathbf{x}, \mathbf{y}, \mathbf{z}] = \text{tcorr}\{\Lambda[\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}], \text{task}[\mathbf{t}]\} = \frac{\text{cov}\{\Lambda[\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}], \text{task}[\mathbf{t}]\}}{\text{std}\_{\text{t}}\{\Lambda[\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}]\} \text{std}\{\text{task}[\mathbf{t}]\}} \tag{14}$$
  $\text{with} \quad \Lambda = \{\text{"A", "P", "J", "J"}^{\text{true}}\}$ 

where std*<sup>t</sup>* 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.
