Open access peer-reviewed chapter

BOLD fMRI Simulation

Written By

Zikuan Chen and Vince Calhoun

Submitted: January 18th, 2016 Reviewed: March 24th, 2016 Published: August 24th, 2016

DOI: 10.5772/63313

Chapter metrics overview

1,695 Chapter Downloads

View Full Metrics


Background: Brain functional magnetic resonance imaging (fMRI) is sensitive to changes in blood oxygenation level dependent (BOLD) brain magnetic states. The fMRI scanner produces a complex-valued image, but the calculation of the original BOLD magnetic source is not a mathematically tractable problem. We conduct numeric simulations to understand the BOLD fMRI model.


  • magnetic resonance imaging (MRI)
  • blood oxygenation level dependence (BOLD)
  • magnetic susceptibility source
  • dipole effect
  • voxelization
  • complex-valued magnetic resonance signal (image)
  • intravoxel dephasing signal
  • multivoxel image
  • BOLD fMRI simulation
  • task correlation

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 [14]. 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 [79] 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.

Figure 1.

A BOLD fMRI model consists of three stages. The stage of “Source Magnetism” provides a dynamic magnetic susceptibility source for the stage of “MRI Technology”. The MRI detection produces a 4D complex-valued fMRI dataset, which are used for functional imaging and mapping by “Statistical Image Analysis”.

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

χ(r,t)=χ0(r)+Δχ(r,t)+noisewith  χ0(r)=χstatic(r)                                                 (parenchymal tissue)Δχ(r,t)=Hctχdo(1Y(t))NAB(r)V(r,t)    (BOLD activity)E1

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

V(r,t)={1,    rvessel0,   otherwise  s.t.  bfrac(t)=(x,y,z)V(r,t)size(V)E2

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

Figure 2.

Illustrations of VOI vasculature and BOLD Δχ source. The VOI is filled with (a1) random vessels (cylindrical segments) and (b1) spheric beads. The NAB-modulated Δχ distributions are shown in (a2) and (b2), respectively, with a y0-slice. It is noted Δχ may assume positive and negative values in local NAB regions.

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 y0-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 μm3, is used to represent a small VOI of 4.1 × 4.1 × 4.1 mm3. 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 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]

b(x,y,z)=B0χ(x,y,z)hdipole(x,y,z) with  hdipole(x,y,z)=14π3z2(x2+y2+z2)3/2(x2+y2+z2)5/2  E3

where * denotes spatial convolution, and hdipole 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]

FT{b(x,y,z}=Hdipole(kx,ky,kz)FT{χ(x,y,z)}B0with Hdipole(kx,ky,kz)FT{hdipole(x,y,z)}=13kz2kx2+ky2+kz2  E4

where (kx,ky,kz) denotes the coordinates in the Fourier domain. The fieldmaps induced by the Δχ distribution in a main field B0 are illustrated in Figure 3 (displayed with a y0-slice), which shows a conspicuous dipole effect in a manifestation of bipolar-valued quadruple lobes around vessels (with an orientation dependence [19]).

Figure 3.

The fieldmaps calculated from the Δχ distributions in Figure 2(a2, b2). It is noted that the Δχ-induced fieldmap 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 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]

C[x,y,z;X]=1|Ω|(x',y'.z')Ω(x,y,z)eiγb(x',y',z')TEwith  X={'TE','B0','|Ω|','vesselsize',}E6

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 (B0), 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

C[x,y,z;X]=1Nn=1N<|Ω|eiγb(xn,yn,zn)TE for (xn,yn,zn)Ω(x,y,z)E7

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

{A[x,y,z;X]=1|C[x,y,z;X]|       (magnitude loss)  P[x,y,z;X]=C[x,y,z;X]           (phase accrual)E8

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

ΩIV(x,y,z)={1    (x,y,z)vessel0      otherwise     ΩEV(x,y,z)={1      (x,y,z)vessel0            otherwise   =1ΩIV(x,y,z)with    Ω=ΩIVΩEV and  |Ω|=|ΩIV|+|ΩEV|E9

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

CIV[x,y,z;X]=1|ΩIV|(x',y'.z')ΩIV(x,y,z)eiγb(x',y',z')TECEV[x,y,z;X]=1|ΩEV|(x',y'.z')ΩEV(x,y,z)eiγb(x',y',z')TEwith    Ω=ΩIVΩEVE10

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.

Figure 5.

Illustration of extravascular (EV) and intravascular (IV) space partition in a voxel for intravoxel spin dephasing 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.

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/r3) 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

C[x,y,z;X]=bfracCIV[x,y,z;X]+(1bfrac)CEV[x,y,z;X]with    bfrac=|ΩIV||Ω|E11

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]

r(t)=x(t)x^+y(t)y^+z(t)z^x(t+δt)=Gauss(x(t),σd)y(t+δt)=Gauss(y(t),σd)z(t+δt)=Gauss(z(t),σd)with   σd=2DδtD={1.5×105cm2/s    (x,y,z)vessel0.75×105cm2/s  (x,y,z)vessel E12

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.

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 timeseries 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

Λtcorr[x,y,z]=tcorr(Λ[x,y,z,t],task[t])cov(Λ[x,y,z,t],task[t])stdt(Λ[x,y,z,t])std(task[t])with   Λ={'A', 'P', 'χtrue', 'χrecon'}E14

where stdt 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 (TE=[0, 60] ms) and for a range of field strength (B0=[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)).

Figure 7.

Intravascular (IV) and extravascular (EV) voxel signal simulations. The IV signal evolves drastically for a long TE and its contribution to the full voxel signal is relatively small.

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

Figure 8.

Multiresolution complex-valued voxel signals due to voxel subdivision. As the voxel size is dyadically reduced, the smaller voxels contain less vessels, and the voxel signal may become turbulent at vessel boundary (Adapted from [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 TE = [0, 60] ms with different field strengths (in terms of ΔχB0 = [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.

Figure 9.

Effects of diffusion and field strength on voxel signal magnitude and phase. It is seen that the diffusion has more effects on magnitude signal than on phase signal and that the diffusion effect decreases as the field strengths increases (Adapted from [20]).

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.

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.

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 B0//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

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 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 y0-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,tON]; (c) an OFF state χ[r,tOFF]; 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, B0 = 3 T, VOI matrix = 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.

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.

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

Numerical simulations of fmap extractions from 4D BOLD fMRI datasets in the presence of additive Gaussian 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 Atcorr for BOLD magnitude fmap from the 4D magnitude image dataspace, and a 3D Ptcorr for BOLD phase fmap from the 4D phase image dataspace. By repeating the 4D BOLD fMRI simulations with different noise levels, we show that Atcorr or Ptcorr is sensitive to the additive Gaussian noise. In Figure 15 are showed the Atcorr and Ptcorr (displayed with a y0-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 (Atcorr 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 Ptcorr) 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 (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-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, [79], 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 [911]. 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 B0, 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):

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

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

  3. As voxel size decreases, the voxel signals evolve more drastically and turbulently inside and around the large vessels.

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

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



1D: one dimensional; 2D: two dimensional; 3D: three dimensional (spatial); 4D: four dimensional (spatiotemporal); BOLD: blood oxygenation level dependent; MR: magnetic resonance; MRI: magnetic resonance imaging; fMRI: functional MRI; FFT: fast Fourier transform; NAB: neuroactive blob; VOI: volume of interest; IV: intravascular; EV: extravascular; gridel: grid element; bfrac: blood volume fraction; fmap: functional map; corr: correlation (coefficient); tcorr: temporal correlation or task correlation.


  1. 1. J. L. Boxerman, L. M. Hamberg, B. R. Rosen et al., “MR contrast due to intravascular magnetic susceptibility perturbations,” Magn Reson Med, vol. 34, no. 4, pp. 555-66, Oct, 1995.
  2. 2. S. Ogawa, R. S. Menon, D. W. Tank et al., “Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model,” Biophys J, vol. 64, no. 3, pp. 803-12, Mar, 1993.
  3. 3. S. Ogawa, D. W. Tank, R. Menon et al., “Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging,” Proc Natl Acad Sci U S A, vol. 89, no. 13, pp. 5951-5, Jul 1, 1992.
  4. 4. O. J. Arthurs and S. Boniface, “How well do we understand the neural origins of the fMRI BOLD signal?” Trends Neurosci, vol. 25, no. 1, pp. 27-31, Jan, 2002.
  5. 5. E. M. Haacke, R. Brown, M. Thompson et al., Magnetic Resonance Imaging Physical Principles and Sequence Design, John Wiley & Sons, Inc, 1999.
  6. 6. Z. Chen and V. Calhoun, “Nonlinear magnitude and linear phase behaviors of T2* imaging: theoretical approximation and Monte Carlo simulation,” Magn Reson Imaging, vol. 33, no. 4, pp. 390-400, May, 2015.
  7. 7. J. L. Boxerman, P. A. Bandettini, K. K. Kwong et al., “The intravascular contribution to fMRI signal change: Monte Carlo modeling and diffusion-weighted studies in vivo,” Magn Reson Med, vol. 34, no. 1, pp. 4-10, Jul, 1995.
  8. 8. J. Martindale, A. J. Kennerley, D. Johnston et al., “Theory and generalization of Monte Carlo models of the BOLD signal source,” Magn Reson Med, vol. 59, no. 3, pp. 607-18, Mar, 2008.
  9. 9. A. P. Pathak, B. D. Ward, and K. M. Schmainda, “A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: the finite perturber method,” Neuroimage, vol. 40, no. 3, pp. 1130-43, Apr 15, 2008.
  10. 10. Z. Chen, and V. Calhoun, “Volumetric BOLD fMRI simulation: from neurovascular coupling to multivoxel imaging,” BMC medical imaging, vol. 12, pp. 8, 2012.
  11. 11. Z. Chen, and V. Calhoun, “Understanding the morphological mismatch between magnetic susceptibility source and T2* image,” Magnetic Resonance Insights, vol. 6, pp. 65-81, 2013.
  12. 12. Z. Chen, A. Caprihan, and V. Calhoun, “Effect of surrounding vasculature on intravoxel BOLD signal,” Medical physics, vol. 37, no. 4, pp. 1778-87, Apr, 2010.
  13. 13. W. M. Spees, D. A. Yablonskiy, M. C. Oswood et al., “Water proton MR properties of human blood at 1.5 Tesla: magnetic susceptibility, T(1), T(2), T*(2), and non-Lorentzian signal behavior,” Magn Reson Med, vol. 45, no. 4, pp. 533-42, Apr, 2001.
  14. 14. Z. Chen, and V. D. Calhoun, “Magnitude and phase behavior of multiresolution BOLD signal,” Concepts Magn Reson, vol. 37B, no. 3, pp. 129-35, 2010.
  15. 15. Z. Chen, and V. Calhoun, “A computational multiresolution BOLD fMRI model,” IEEE Trans. BioMed. Eng, vol. 58, no. 10, pp. 2995-9, 2011.
  16. 16. D. A. Yablonskiy, “Quantitation of intrinsic magnetic susceptibility-related effects in a tissue matrix. Phantom study,” Magn Reson Med, vol. 39, no. 3, pp. 417-28, Mar, 1998.
  17. 17. J. R. Reitz, F. J. Milford, and R. W. Christy, Foundations of Electromagnetic Theory, New York: Addison-Wisley, 1993.
  18. 18. J. P. Marques and R. Bowtell, “Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility,” Concepts Magn. Reson, vol. B 25, pp. 65-78, 2005.
  19. 19. Z. Chen and V. Calhoun, “Effect of object orientation angle on T2* image and reconstructed magnetic susceptibility: numerical simulations,” Magn. Reson. Insight, vol. 6, pp. 23-31, 2013.
  20. 20. Z. Chen and V. Calhoun, “Computed diffusion contribution in the complex BOLD fMRI signal,” Conc. Magn. Reson. Part A, vol. 40A, no. 3, pp. 128-145, 2012.
  21. 21. Z. Chen and V. Calhoun, “Computed inverse resonance imaging for magnetic susceptibility map reconstruction,” J. Comp. Assist. Tomo, vol. 36, no. 2, pp. 265-74, Mar, 2012.
  22. 22. Z. Chen and V. Calhoun, “Reconstructing brain magnetic susceptibility distributions from T2* phase images by TV-reguarlized 2-subproblem split Bregman iterations,” Rep. Med. Imag, vol. 7, pp. 41-53, 2014.
  23. 23. Z. Chen, Z. Chen, and V. Calhoun, “Blood oxygenation level-dependent functional MRI signal turbulence caused by ultrahigh spatial resolution: numerical simulation and theoretical explanation,” NMR Biomed, vol. 26, no. 3, pp. 248-64, Mar, 2013.
  24. 24. P. Latta, M. L. Gruwel, V. Jellus et al., “Bloch simulations with intra-voxel spin dephasing,” J. Magn. Reson, vol. 203, no. 1, pp. 44-51, Mar, 2010.

Written By

Zikuan Chen and Vince Calhoun

Submitted: January 18th, 2016 Reviewed: March 24th, 2016 Published: August 24th, 2016