**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 transfor‐ mations (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.
