*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

**Figure 10.** Illustration of VOI configuration and voxelization. A VOI is represented by a large matrix for subvoxel structure representation. The VOI voxelization produces a small multivoxel matrix, depending on the voxel size.

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 voxeli‐ zation 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.

**Figure 11.** Illustration of volumetric BOLD fMRI simulation, displayed with a *z*-slice (with *B*0//*z*-axis). (a) 3D Δχ source (in a matrix of 64 × 64 × 64 voxels resulting from a VOI of 2048 × 2048 × 2048 gridels); (b) the Δχ-induced fieldmap; (c, d) the magnitude and phase images with static spins (at *TE*= 30 ms); and (e, f) the magnitude and phase images with diffusion spins.

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

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 voxeli‐ zation 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

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

**Figure 11.** Illustration of volumetric BOLD fMRI simulation, displayed with a *z*-slice (with *B*0//*z*-axis). (a) 3D Δχ source (in a matrix of 64 × 64 × 64 voxels resulting from a VOI of 2048 × 2048 × 2048 gridels); (b) the Δχ-induced fieldmap; (c, d) the magnitude and phase images with static spins (at *TE*= 30 ms); and (e, f) the magnitude and phase images with

size is comprised of more spatial smoothing.

18 Numerical Simulation - From Brain Imaging to Turbulent Flows

*3.2.2. Multivoxel image calculation*

insignificant on image formation.

diffusion spins.

We performed volumetric BOLD fMRI simulations for a span of echo times (*TE*= [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.

**Figure 12.** Spatial correlation measurements (a) between χ source and magnitude image and (b) between χ-induced fieldmap and phase image, for static intravoxel dephasing and diffusive intravoxel dephasing. Note the small display windows for corr values in the range of [–1,1] (Adapted from [11]).

#### **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 NABmodulated 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).

**Figure 13.** Numerical representation of a local BOLD activity in terms of 4D χ(**r**,*t*). (a) A Gaussian-shaped NAB and a ball-shaped NAB in VOI; (b) an ON state χ[**r**,*t*ON]; (c) an OFF state χ[**r**,*t*OFF]; and (d) the bead-laden structure in a voxel.

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 (*TE*= 30 ms, *B*0 = 3 T, VOI matrix =

**Figure 14.** Numerical simulations of 4D BOLD fMRI data acquisition. (a1, a2) BOLD magnitude images captured at an ON and OFF state and (a3) the magnitude signal timecourses at two voxels (marked by x and *o* (a1, a2), extracted from the 4D magnitude dataset *A*[**r**,*t*]). (b1, b2, b3) for the BOLD phase images in *P*[**r**,*t*] and the voxel phase timecourses.

64 × 64 × 64 voxels, 1 voxel = 32 × 32 × 32 gridels, 1 gridel = 2 × 2 × 2 μm3 ). **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).

**Figure 13.** Numerical representation of a local BOLD activity in terms of 4D χ(**r**,*t*). (a) A Gaussian-shaped NAB and a ball-shaped NAB in VOI; (b) an ON state χ[**r**,*t*ON]; (c) an OFF state χ[**r**,*t*OFF]; and (d) the bead-laden structure in a voxel.

20 Numerical Simulation - From Brain Imaging to Turbulent Flows

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 (*TE*= 30 ms, *B*0 = 3 T, VOI matrix =

**Figure 14.** Numerical simulations of 4D BOLD fMRI data acquisition. (a1, a2) BOLD magnitude images captured at an ON and OFF state and (a3) the magnitude signal timecourses at two voxels (marked by x and *o* (a1, a2), extracted from the 4D magnitude dataset *A*[**r**,*t*]). (b1, b2, b3) for the BOLD phase images in *P*[**r**,*t*] and the voxel phase timecourses.

**Figure 15.** Numerical simulations of fmap extractions from 4D BOLD fMRI datasets in the presence of additive Gaussi‐ an noises at different noise level = {0.001,0.01,0.05,0.1}. (a) Magnitude fmap and (b) phase fmap.

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 spatiotem‐ 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 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.
