PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Concepts Magn Reson Part B Magn Reson Eng. Author manuscript; available in PMC 2010 September 29.
Published in final edited form as:
PMCID: PMC2947032
NIHMSID: NIHMS194758

Magnitude and phase behavior of multiresolution BOLD signal

Abstract

High spatial resolution fMRI provides a more precise estimate of brain activity than low resolution fMRI. The magnitude and phase parts of the BOLD signals are impacted differently by changes in the scan resolution. In this paper, we report on a numerical simulation to show the impact of spatial resolution upon the complex-valued BOLD signal in terms of magnitude and phase variation. We generate realistic capillary networks in cortex voxels, calculate the BOLD-induced magnetic field disturbance and the complex BOLD signals for the voxel and its subvoxels, and thereby characterize the magnitude and phase behaviors across multiple grid resolutions. Our results show that: 1) at higher spatial resolution there is greater spatial variation in the phase of the BOLD signal as compared to its magnitude; 2) the spatial variation of the phase signal monotonically increases with respect to spatial resolution while for the magnitude the spatial variation may reach a maximum at some resolution level; 3) voxels containing large capillaries have higher phase spatial variation than those with smaller capillaries; 4) the amplitude spatial variation at a resolution level increases with respect to relaxation time whereas the phase variation is generally unaffected.

Keywords: Computational fMRI, BOLD simulation, high spatial resolution, voxel subdivision, BOLD signal, magnitude, phase, contrast, gridel

I. Introduction

For accurate spatial depiction of anatomical structure and physiological functions, there is an ongoing push towards higher spatial resolution functional MRI (fMRI) 1-7 Higher resolution fMRI will improve the brain activity mapping in terms of spatial localization, topographic characterization, and topological connections of neuron activities. The blood-oxygenation-level-dependent (BOLD) signal is attributed to spin procession of protons in BOLD-induced inhomogeneous magnetic field (called BOLD field henceforth). Since the spin precession of a single proton is technically undetectable with the current MRI technology, what we observe in fMRI is the collective behavior of a tremendous number of protons in individual voxels 8. The BOLD signal resulting from intravoxel spin dephasing is represented by a complex number consisting of amplitude and phase parts. However, only its amplitude part is widely utilized; the phase part is not found in use so far. Due to the detection of collective spin precession, a voxel BOLD signal conceals the individual behaviors of protons within the voxel. By spatial subdivision, a voxel is partitioned into subvoxels; correspondingly, the voxel signal is decomposed into subvoxel signals thus revealing additional spatial information. In this paper, we show BOLD contrast regeneration by recursively decomposing voxels and BOLD signals.

Currently, the phase part of the fMRI signals is not typically used. One reason for this may be due to unrestricted phase accumulation and phase wrapping associated with the complex voxel signal formation9-11. Another reason may be due to the fact that the amplitude pattern (assuming a nonnegative intensity distribution) is intuitively appreciable because it reproduces the realistic scene, whereas the phase pattern (assuming both positive and negative distributions) manifests itself in an encrypted form (in terms of spatial or temporal shifts) and consequently is not straightforwardly perceivable.12,13 Perhaps the primary reason for phase is ignored is due to the lack of fMRI phase mechanisms, as evidenced from the ongoing research6,9-11,14-23. Toward the goal of making use of BOLD signal phase, in this paper, we attempt to understand the BOLD phase behavior at different spatial resolutions by numerical simulation.

II. Methodology

II.1 Voxel partition and gridels

An fMRI image reflects the brain activity with a finite spatial resolution (~mm). Brain anatomy shows that the brain cortex consists of capillary vasculature (radius [set membership] [2, 10] micron), which is responsible for BOLD activity (typically causing magnetic field perturbation < 1ppm: parts per million). We can simulate a cortex voxel by filling it with a random capillary vasculature under the condition of a fixed blood volume fraction. Let V0 denote a cortex voxel of size D0×D0×D0 (D0~mm), which can be expressed by a binary volume

V0(x,y,z)={1,(x,y,z)vessel0,otherwisexD02,yD02,zD02,
[1]

As spatial resolution increases, the voxel size decreases. The multigrid voxel subdivision scheme is depicted in Fig. 1(a), which is an octadic spatial partition by a 2×2×2 grid, producing 8 subvoxels that can be labeled by indices in an octet set K8={1,2,3,4,5,6,7,8}. The recursive voxel partition produces a full octadic tree, octree for short 24, as diagramed in Fig. 1(b). In an octree, a node at level j can be uniquely designed by a j-tuple (1k,2k,…,jk), such as the exemplary illustration of an node at (13,28,37) in Fig. 1(b). In the appendix, we provide the notation and formulas for describing an octree.

Fig 1
Recursive octadic voxel partition, a gridel is a subvoxel (or a bin in the voxel), which can be used as a mask. (b) Graphic octree. A node at resolution level j is uniquely labeled by a j-tuple (1k,2k,…,jk), which traces its path to the root. ...

In principle, a digital image can be reconstructed into (or represented by) a digital array with an arbitrary grid resolution 25. We coin the word “gridel” to designate a grid element that is generated by rectilinear gridding partition. The application of an octadic partition to a gridel produces an octet of children gridels in Fig. 1(a), in which V0 is considered a gridel that has no parents as well. The gridels generated from multigrid partition can be systematically organized by an octree as shown in Fig. 1(b). By assigning each octree node a property and using the j-tuple label, we can construct a multigrid volume (or mask) octree Vj(1k,2k,…,jk), a multigrid BOLD field octree Bj(1k,2k,…,jk), and a multigrid BOLD signal octree Sj(1k,2k,…,jk), and so forth. We also call a multigrid octree a gridel octree because the multigrid decomposition produces a collection of gridels that are used to construct an octree. In particular, the gridel at level j takes on a cubic bin in size of dj× dj× dj, with dj=D0/2j, which can be considered as a local mask in voxel subdivision. As a result, the voxel representation and manipulation for multiresolution BOLD simulation involve three important masks: Vε, Vj and V0, which are governed by a spatial containment relationship: Vε [subset, dbl equals] Vj [subset, dbl equals] V0. Specifically, Vε is the building block for both voxel V0 and gridels Vj; V0 is initially supported or defined by as many as Nε 3 building blocks (fixed number); a gridel Vj consists of a number of Nj 3 building blocks, usually, 1 < Nj [double less-than sign] 2Nε for 0 < j [double less-than sign] Nε; in extreme cases, a gridel Vj could be as large as V0 at j=0 or as tiny as Vε at j=Nε. Based on a finest grid array representation (in size of Nε×Nε×Nε) of a voxel V0, our multiresolution voxel and signal decomposition strategy is similar to the “coarse-to-fine” scheme used in wavelet transform26.

II.2 Multigrid BOLD field partition

Given the vessel geometry for a cortex voxel, we assume that blood flow is confined in vessels, and the brain activity changes the oxygenation level in the blood stream in the cortex region. This BOLD mechanism allows us to express the 3D blood susceptibility variation pattern by 12

Δχ(x,y,z)=ΔχbV0(x,y,z)withΔχb=(1Y)ΔχdoHct
[2]

where Hct stands for blood hematocrit, Y denotes oxygenation saturation level, Δχdo the susceptibility difference between oxygenated and deoxygenated blood. Eq. [2] serves as the connection between our numerical simulation and biology (vessel network) and physiology (blood susceptibility variation).

From the susceptibility distribution, we can calculate the BOLD field distribution 27-29

ΔB(x,y,z)=IFT{(kz2kx2+ky2+kz213)FT{V0(x,y,z)}}ΔχbB0
[3]

where FT and IFT stands for Fourier transform and inverse Fourier transform, respectively, and (kx, ky, kz) represents the Cartesian coordinates in Fourier domain. In fact, the z-component of ΔB(x,y,z) is the blood-susceptibility-induced magnetic field, which is responsible for BOLD signal.

In conformation with the support grid of the voxel V0, the spatial distribution of ΔB over V0 can also be discretized on the Nε×Nε×Nε array by using the mask bins in Eq. [A2]. Corresponding to octreeJ(V0) in Eq. [A4], we can construct a BOLD field tree by calculating the gridel-averaged BOLD field by

ΔBj(k1,k2,,kj)=ixiyizΔB(ix,iy,iz)Vj(ix,iy,iz;k1,k2,,kj)Nj3
[4]

where (ix, iy, iz) denotes the discrete locations inside the voxel (see Eq. [A1]). There are 8j entries for ΔBj at level j, we can perform statistics in terms of {‘standard deviation (std)’, ‘maximum (max)’, ‘minimum (min)’} to describe the local BOLD field inhomogeneity. Meanwhile, for visualization purposes, we can arrange the 8j entries into an Mx×My matrix under the condition of MxMy=8j. Specifically, the matrix arrangements for 4-level decomposed entries are as follows: 2×4 (=8), 8×8 (=82), 32×16 (=83), and 64×64 (=84).

II.3 Multigrid BOLD signals

Due to the indivisibility of Vε, the magnetic field in Vε is assumed to be locally uniform (with a single value in numerical simulation). In MRI physics, a single spin in the magnetic field ΔB can be expressed by exp(-iγΔBt), where γ is the gyromagnetic ratio and t the relaxation time 12. Albeit a small macroscopic dimension (~ micron), Vε may contain a tremendous number of protons (counted in units of mole or Avagadro’s number), which are indistinguishable to signal detection by the current MRI technology. Thus we model the collective precession behavior in Vε at (ix,iy,iz) by 23

Sε(ix,iy,iz,t)=exp(iγΔB(ix,iy,iz)t))
[5]

Although there are numerous protons in Vε, we calculate the signal for Vε in Eq.[5] by a single spin precession in magnetic field ΔB(ix,iy,iz). The agglomerative spin precessions at a position in BOLD field has been described by a concept of “spin packet” 30, for which we will implement it digitally by an illustration in Fig. 2. The “spin packet” model allows us to calculate the BOLD signal resulting from numerous protons in gridel Vj by manipulating a manageable number of Nj 3 spin packets. On the other hand, the “spin packet” model sets the limit of simulation accuracy due to the indivisibility in a spin packet. Exposed to inhomogeneous BOLD field, the spin packets distributed in a gridel Vj precess out of phase, which is described as precession dephasing that is responsible for gridel BOLD signal formation.

Fig 2
Illustration of the “spin packet” model and “gridel”. The collective spin precession of numerous protons in Vε is modeled as a single spin packet; a gridel (marked in red box) assumes a bin within a voxel and it ...

The “spin packet” model enables the computation implementation (by Eq. [5]) for the aggregation of numerous protons that render spin precessions in the local magnetic field over the build block Vε. It is a practical approximate model due to its finite size. This model can be justified when Vε is substantially tiny such that the local magnetic field over Vε is uniform and the spins in Vε render in-phase procession. In theory, when ε approaches infinitesimal, ε→0, the “spin packet” model will eventually reduce to the single proton model; nevertheless, the astronomical number of protons are beyond the enumeration by the state-of-the-art computer. In practice, ε[double less-than sign]D0, the current MRI technology limits the spatial resolution (V0[dbl greater-than sign] Vε), therefore, we don’t need to simulate the spatial resolutions of BOLD signal far beyond the experimental feasibility. Suppose a voxel is in size of 320×320×320 micron3, for example, the MRI resolution is 20 micron, then 4-level decomposition suffices (see our simulation later: j={1,2,3,4} and Nε=512). So we only need to render a few levels (j [double less-than sign] Nε) for the coarse-to-fine multigrid decomposition. As discussed above, the number of gridel nodes in an octree proliferates exponentially (by 8j, see Eq. [A2]). Let Sj(t;1k,2k,…,jk) denote the BOLD signal generated from all spin packets within a gridel at (1k,2k,…,jk), then it can be calculated by

Sj(t;k1,k2,,kj)=ixiyizSε(ix,iy,iz,t)Vj(ix,iy,iz;k1,k2,,kj)Nj3
[6]

where Vj serves as a spatial mask (in Eq. [A2]) that only counts the spin packets within the gridel in question.

Usually, we need to examine the amplitude and phase of the BOLD signals at level j by performing a nonlinear operation on Sj(t), that is,

Aj(t;k1,k2,,kj)=Sj(t;k1,k2,,kj)Φj(t;k1,k2,,kj)=argSj(t;k1,k2,,kj)withSj(t;k1,k2,,kj)=Aj(t;k1,k2,,kj)exp(iΦj(t;k1,k2,,kj))
[7]

where the complex BOLD signal is expressed by amplitude Aj(t) and phase Φj(t).

II. 4. Characterization of multiresolution amplitude and phase behavior

In appendix B, we depict the multiresolution BOLD signal tree and its properties. Based on the octree notations, we can characterize the amplitude variation and phase variation behavior not only at a resolution level but also across multiple resolution levels.

The multigrid BOLD signals possess global and local rebinning invariance (see Eqs. [B3] and [B4]) if we do not interfere with them nonlinearly. In practice, we must get the amplitude phase parts separately from a complex-valued BOLD signal by Eq. [7]. It is the nonlinear operation that brings out the inequalities in Eqs.[B5] and [B6], which are responsible for regenerating the signal differences across gridels (at the same level j). The number of spatially decomposed signal is dependent upon the resolution level, and the dispersion among these entries at a level depends on the relaxation time as well. The dispersion behavior can be statistically characterized by the coefficient of variation (defined by a ratio of std/mean) 31. With the reference to non-decomposed voxel signal, S0(t)=A0(t)exp(-iΦ0(t)), we characterize the amplitude and phase dispersions respectively by their coefficients of variations (denoted by cv), that is

cvAj(t)=std{Aj(t;k1,k2,,kj)k1,k2,,kjK8}A0(t)
[8]

cvPj(t)=std{Φj(t;k1,k2,,kj)k1,k2,,kjK8}Φ0(t)
[9]

It is noted that the dimensionless ratios permits the comparison of cvA and cvP with the same measure (dimensionless scale).

At each level j, there are as many as 8j entries associated with a parameter such as the gridel signal. For visualization purposes, we can arrange them into an Mx×My matrix under the layout of MxMy= 8j, in the same way done for ΔBj display. When the amplitude entries are rearranged into an Mx×My matrix, we can examine its non-uniformity statistically in terms of {‘std’, ‘max’, ‘min’} in the sense of image analysis. In order to include the comparison with the non-decomposed voxel signal S0(t) in the matrix-formatted dispersion characterization, we define the amplitude contrast (conA) and phase contrast (conP) respectively by

conAj(t;k1,k2,,kj)=Aj(t;k1,k2,,kj)A0(t)A0(t)
[10]

conPj(t;k1,k2,,kj)=Φj(t;k1,k2,,kj)Φ0(t)Φ0(t)
[11]

Statistically, “std” can represent the statistical variation, and “max” the extremity (or outlier) distribution. In order to characterize the unpredictable variation of the entries generated at multiresolution levels with both “std” and “max”, we define a dimensionless parameter, in terms of “wildness”32, by

wildAj(t)=max{Aj(t;k1,k2,,kj)A¯j(t)k1,k2,,kjK8}std{Aj(t;k1,k2,,kj)k1,k2,,kjK8}A¯j(t)=k1K8k2K8kjK8Aj(t;k1,k2,,kj)8j
[12]

wildPj(t)=max{Φj(t;k1,k2,,kj)Φ¯j(t)k1,k2,,kjK8}std{Φj(t;k1,k2,,kj)k1,k2,,kjK8}Φ¯j(t)=k1K8k2K8kjK8Φj(t;k1,k2,,kj)8j
[13]

Again, the dimensionlessness of wildA and wildP allows us to compare the outlier excursions of amplitude and phase with the same measure.

In a summary, we calculate the time course signals (with a relaxation time t [set membership] [0,120]ms) for all the gridels, thereby to reveal the amplitude and phase variation behaviors across multiple resolution levels by {cvAj(t), cvPj(t)} and {wildAj(t), wildPj(t)} with j={1,2,3,4}, and the behaviors at a resolution level by {conAj(t; 1k,2k,…,jk), conPj(t; 1k,2k,…,jk)} with a fixed j. All the parameters reduce to zeros for the non-decomposed voxel, that is, cvA0=0, cvP0=0, conA0=0, conP0=0, wildA0=0, wildP0=0 because of the singleton of S0(t). For most cases, the following inequalities are observed:

cvPj(t)>cvAj(t)wildPj(t)>wildAj(t)std{conPj(t;k1,k2,,kj)k1K8,,kjK8}>std{conAj(t;k1,k2,,kj)k1K8,,,jkK8}
[14]

cvPj+1(t)>cvPj(t)wildPj+1(t)>wildPj(t)std{conPj+1(t;k1,k2,,kj)k1K8,,kjK8}>std{conPj(t;k1,k2,,kj)k1K8,,,jkK8}
[15]

Specifically, the inequalities in Eq.[14] depict the amplitude and phase variations compared at the same resolution level, and those relations in Eq.[15] reveal the across-resolution-level behavior for phase part, which are not necessarily valid for magnitude part. Although, these relations are beyond analytic proof, which are of practical value for noise characterization or high contrast fMRI pursuit. The goal of this paper is to quantitatively characterize these inequalities by numeric simulation, as reported in next section.

III. Numeric simulation

The BOLD fMRI simulation requires a collection of parameters and their settings 21,23,33,34. Table 1 lists the parameter settings used in our simulation.

Table 1
Parameters for numerical simulation

We generated cortex voxels (D0=320 micron) in a support grid of 512×512×512 with capillaries (radius = {2.5, 6, 9} micron) under the condition of 2% blood volume fraction. Three voxel geometries are shown in Fig. 3.

Fig. 3
Cortex voxels filled with three capillaries: radius={2.5, 6, 9}micron, with the same blood volume fraction of 0.02.

For each voxel V0, we performed 4-level multigrid decomposition to generate an octree, octree4(V0), as a realization of Eq.[A4]. Table 2 collects the octree properties: each gridel at level j assumes a size of dj×dj×dj, and contains a number of Nj 3 spin packets; the total number of gridels at level j is 8j. The level j=Nε is included in the rightmost column as the simulation limit (ε = 0.625 micron), which could not be reached by the current MRI technology.

Table 2
Properties and values of an octree.

Given a voxel geometry, we simulated the BOLD activity with the parameter setting in Table 1, calculated the BOLD-induced magnetic field by Eq.[3], and performed multigrid decomposition on the BOLD field with the grid masks in octree4(V0). Figs. 4 through through1111 only demonstrate the results done for the cortex voxel of 6-micron capillary. The simulation results for vessel networks of three capillaries (radius={2.5, 6, 9}micron) are collected in tables and are presented in Figs 12 through through1414.

Fig. 4
Demonstration of multigrid BOLD magnetic field partition at a slice ΔBj(ix,iy,iz=20) for : partitions of (a) 2×2 for j=1, (b) 4×4 for j=2, (c) 8×8 for j=3. The support grid for a gridel at level j=3 is magnified in (d). ...
Fig. 11
Matrix representation of phase contrasts (conPj defined in Eq.[11]) of multigrid BOLD signals. The entry value range is indicated in the grayscale bar (unit: dimensionless).
Fig. 12
Amplitude and phase variation behaviors due to multiresolution BOLD signal refinement (t=30ms). The upper row shows the variation behaviors in terms of standard deviation (std), the middle row in coefficient of variation (cvA and cvP in Eqs. [8] and ...
Fig. 14
Effect of vessel capillary size on magnitude and phase variations during multiresolution BOLD signal refinement (t=30ms). The resolution levels are indicated by the legend in the left upper panel. Three capillary radiuses ({2.5, 6, 9} micron) are provided ...

Fig. 4 shows a slice of the BOLD field pattern under multigrid partition. In particular, Fig. 4(d) shows a gridel at level j=3, as will as a blowup revealing its support grid mesh. Fig. 5 displays the gridel-averaged multigrid magnetic field patterns in image format (by arranging 8j entries into an Mx×My matrix). It is seen the inhomogeneity of the BOLD field emerges and enhances as the grid resolution refines. The numeric characteristics {‘mean’, ‘std’, ‘max’, ‘min’} for the multigrid BOLD field inhomogeneity for three voxels are collected in Table 3. It shows that the extreme values {‘max’, ‘min’} differ drastically (in magnitude) during multiresolution refinement, which are eventually bounded by the support grid level (j=Nε). This BOLD signal refinement behavior implies that there exist considerable cancellation of BOLD field superimposition within gridels.

Fig. 5
Matrix presentation of multigrid BOLD magnetic field (spatial subdivision of blood-susceptibility-induced magnetic field inhomogeneity over voxel V0), calculated by gridel average by Eq. [4] for: (a) j=1, (b) j=2, (c) j=3, (d) j=4. The unit of the common ...
Table 3
Inhomogeneity characteristics of multigrid BOLD fields (unit:10-3 ppmT). The means are j-invariant, and other statistics are j-dependent.

Before launching multigrid decomposition, we generated the BOLD signal for V0, that is S0(t)=A0(t)exp(-iΦ0(t)). Fig. 6 shows the amplitude decay (A0(t)) and phase angle accumulation (Φ0(t)) for relaxation time t [set membership] [0, 120]ms. It is seen that the phase accumulates linearly for the overall voxel V0.

Fig. 6
The amplitude and phase of the voxel BOLD signal: S0(t)=A0(t)exp(-iΦ0(t)) (defined in Eq. [7] with j=0) which is generated from the plenary of spin packets (in a number of 1.34×108=512×512×512).

We started to decompose the voxel recursively and calculate the gridel BOLD signals. Based on voxel partition and BOLD field partition, we calculated the gridel BOLD signal Sj(t) by Eq. [6], then extracted amplitude and phase parts by nonlinear operations in Eq. [7]. For 4-level decomposition, the amplitude parts are shown in Fig. 7 and the phase parts in Fig. 8. It is noted that the BOLD signal decay (amplitude) tends to oscillate in Fig. 7(d), and the phase angle accumulates curvedly in Figs. 8(a) and (b) and some entries dance (due to local fast phase accumulation) in Figs 8(c) and (d).

Fig. 7
The amplitude decay behavior of multigrid BOLD signal refinement in terms of Aj(t) (calculated by Eq. [7]) for j={1,2,3,4}, in which A0(t) (in thick broken line) is provided as a reference (calculated by [7] with j=0). The curves are colored for visual ...
Fig. 8
The phase accumulation behavior of multigrid BOLD signal refinement in terms of Φj(t) (calculated by Eq. [7]) for j={1,2,3,4}, in which Φ0(t) (in thick broken line) is provided as a reference (calculated by [7] with j=0). The curves are ...

Based on Eqs. [10] and [11], we calculated conA and conP for the multigrid BOLD signals at each level for a relaxation time range of [0, 120] ms. At relaxation times of t={30, 60, 90}ms, we examined the variations of phase and magnitude across multigrid levels (j={0,1,2,3,4}) in terms of standard deviation with respect to relaxation time and resolution level. The results are shown in Fig. 9, which shows that the amplitude variation at a resolution level increases with respect to the relaxation time, while the phase variation seems not to be affected as the relaxation time lengthens. Fig. 9 also shows that the amplitude variation at a relaxation time tends to be saturated after level j=3, while the phase variation constantly increases with respect to the resolution level. To inspect the contrast regeneration due to multiresolution spatial decomposition, we show conA and conP in matrix format in Figs. 10 and and1111 for 4-level multigrid decomposition and t=30ms.

Fig. 9
Amplitude and phase variation behaviors of multigrid BOLD signal refinement in terms of std{conAj(t)} and std{conPj(t)}, with respect to relation time (t) and multigrid level (j). The inset in lower right panel is the close-up of std{conPj(t)} at j=3. ...
Fig. 10
Matrix representation of amplitude contrasts (conAj defined in Eq.[10]) of multigrid BOLD signals. The entry value range is indicated in the grayscale bar (unit: dimensionless).

The BOLD signal refinement behaviors for three capillary networks are presented in Figs. 12 through through14,14, in terms of {‘std’, ‘cvA’, ‘cvP’, ‘wildA’, ‘wildP’} with respect to three parameters: resolution level j, relaxation time t, and capillary size. Meanwhile, more numerical results for cortex voxels of capillaries of {2.5, 6, 9}micron are collected in Table 4.

Table 4
Numeric simulation results of multigrid BOLD signal refinemen.

Looking at Figs. 9, ,1212 through through14,14, and Table 4, we observe the following behaviors.

  1. Phase variation is much higher than amplitude variation. Specifically, based on the setting of {voxel: 320×320×320 micron3, support grid array: 512×512×512, capillary radius: 6 micron, relaxation time: 60ms, decomposition resolution levels: {1,2,3}} we found an amplitude variation of 0.25 (dimensionless) and a phase angle variation of 0.47 rad, and found that the regenerated phase contrast is about 4 times higher than the amplitude counterpart, and extremely, the phase may excurse more wildly than the amplitude (up to ~ 20 times larger, in terms of max/std).
  2. At the first few low resolution levels (j=1,2), both amplitude and phase variations increase. As the resolution level ever increases (j>3) and as the relaxation time lengthens (t>30ms), the magnitude variation may go down, whereas the phase variation increases consistently and drastically.
  3. For a fixed blood volume fraction (2% in this simulation), large capillaries incur more phase variation, however, this observation is not always true for amplitude variation. In comparison, phase variation is higher than amplitude variation in terms of cvP/cvA and wildP/wildA, especially these ratios increase drastically for large capillaries as resolution refines.

Above all, we should point out that, in principle, the BOLD signal for a gridel should be calculated based on the “spin packet” model in Eq.[6], rather than from gridel-averaged BOLD field in Eq. [4]. The similarity among the matrix displays in Figs. 5 and and1111 suggests an empirical approximation: the BOLD image phase pattern may approximate the gridel-averaged BOLD field. Extremely, Fig. 5 assumes the same pattern as Fig. 11 if Vj=Vε.

IV Discussion

Currently, the spatial resolution fMRI is typically larger than 1 millimeter. For numerical simulation, we can adopt an arbitrary grid resolution (the limit of support grid) provided that the computation workload is feasible. To provide meaningful simulation results, we start with a cortex voxel in size of 320×320×320 micron3 cortex and multigriding it down to gridels in size of 20×20×20 micron3. The support grid in the “spin packet” model is a main factor that determines the simulation accuracy and validity. By using the multigrid masks and adopting the “spin packet” model, we can carry out the numerical calculation of multigrid BOLD signals by manipulating only a handful of spin packets in a gridel. Nevertheless, accurate simulation pursues high-grid-resolution depiction of vasculature or tiny spin packets. Due to limitation on computation facility and time, we represent the vortex voxel in a support grid by a 512×512×512 array. For highlighting the multigrid BOLD signal behavior, we only consider static BOLD inhomogeneous field and stationary spin packet precession. The effect of dynamic BOLD field and spin diffusion will accentuate the intra-gridel dephasing, which deserves further investigation.

The reason for the regeneration of amplitude contrast and phase contrast due to multigrid voxel decomposition lies in the nonlinear operations of extracting the amplitude and phase parts from a complex-valued fMRI signal. If we do not perform nonlinear operations, all the spin packets are just conceptually rebinned (regrouped) by gridels, and the associative addition (for scalars or vectors) is just a trivial linear operation, therefore the mean values are invariant to the multigrid decompositions (see Table 3). We describe this phenomenon in terms of rebinning invariance (see Appendix B) and suggest use this property to check the algorithm implementation. By the way, the nonlinear operations associated with the magnitude and phase contrast regenerations provide an instance for the physical principle: a physical state may be changed upon detection.

Both amplitude and phase in the complex BOLD signal evolve with the relaxation time. In general, the amplitude decays due to spin precession dispersion, and it may damply oscillate for long relaxation time; meanwhile the phase angle accumulates, which manifests a linear accrual for low resolution levels (j=1,2) and develops into sublinear or nonlinear progression for high resolution levels (j=3,4). The spatial resolution refinement brings about more details and irregularities to the distribution of the BOLD field perturbation. An emerging peak in the refined BOLD field map accounts for the fast phase accumulation and phase wrapping, thus for the contrast regenerations of multiresolution magnitudes and phases. Statistically, the phase variation is always higher than the amplitude variation 35, and this behavior is expected to be exaggerated for large vessels at fine resolutions.

Due to bipolar BOLD field distribution (assuming both positive and negative values), the phase angles diverse in positive and negative directions during multiresolution decomposition (see Fig.8). We observed that, even for a linear accumulative phase angle for a voxel signal (see Fig. 6), the spatial multiresolution decomposition could spawn a multiple of phase angles that could excurse about Φ0(t), straying in two opposite directions (positive up or negative down). In most practical case, the BOLD signal S0(t) for a cortex voxel produces a non-negative amplitude part (A0(t)[set membership](0,1]) and non-zero phase part (Φ0(t)≠0 for t>0). It is quite unlikely for a cortex voxel to produce a time series of Φ0(t) that stays at or closes to zero during intravoxel dephasing. Furthermore, the time series of A0(t) and Φ0(t) stay compromisingly and smoothly amidst the dispersions of Aj(t) and Φj(t)), respectively. Therefore, we suggest to use S0(t) as the reference for gridel signal analysis at all resolution levels. For the purposes of comparing amplitude and phase dispersions, we always use the dimensionless measures, such as {‘cvA’, ‘cvP’} and {‘wildA’, ‘wildP’}.

V. Conclusion

As spatial grid resolution refines, local inhomogeneity of the BOLD-induced magnetic field, which was once degenerated due to voxel average at a coarse resolution, now becomes detectable. Examining the amplitude and phase parts of the BOLD signal separately, we observed that phase variation is consistently higher than the amplitude variation as the resolution refines. For particular voxels filled with micron capillaries, with the setting of resolution level j=1 and relaxation time t=30ms, we found that the phase variation is larger than the amplitude variation by a factor of about 4, and this ratio continues to rise as the resolution refines. At a resolution level, the amplitude variation increases with respect to relaxation time, whereas the phase variation is basically unaffected in absence of diffusion as the relaxation time lengthens (<120ms). We also found that the amplitude variation during multiresolution refinement may assume a maximum at a certain resolution level for some vascular geometry. Therefore, it is not recommended to pursue higher resolution levels (e.g., j > 3 for the capillary vasculature used in this study) when working only with the amplitude signal. These findings show that magnitude and phase behave differently with respect to {‘spatial resolution’, ‘relaxation time’, ‘capillary size’}. The fact that the phase variation increases with respect to spatial resolution much more than magnitude variation may explain why the phase part in fMRI is usually considered too noisy to be useful, thus remains unexploited so far. On the other hand, if the phase mechanism is well understood and the factors that cause phase contrast can be controlled, we can make use of it for high-contrast BOLD imaging.

Fig. 13
Comparison of amplitude and phase variations due to multiresolution spatial refinement of voxel BOLD signal. The comparison was done with three voxel vessel geometries: radius= {2.5, 6, 9}micron, and 30ms relaxation time. The upper row shows the comparison ...

Acknowledgments

This research was supported by the NIH (1R01EB006841, 1R01EB005846) and NSF (0612076). The authors thank Dr. Arvind Caprihan for helpful discussion.

Appendix A

Octree notation for multigrid voxel decomposition

Based on recursive octadic subdivision, the multigrid voxel partition can be represented by a regular octave tree (octree), in which each node can be uniquely designated by its pathway from the tree root. We use K8={1,2,3,4,5,6,7,8} to denote the octet, and a j-tuple (1k,2k,…,jk) to label a node at level j. It is helpful to look at the illustration for node (13,28,37) in Fig. 1(b). For numerical simulation, we need to discretize the voxel V0 in Eq. [1] into an array, called the support grid for V0, as large as possible in the sense of simulation rationale. Suppose V0 be a cube of D0×D0×D0 in physical size (e.g., D0 =320 micron), and we discretize it into Nε×Nε×Nε array (Nε[dbl greater-than sign]1) with the grid resolution ε by

{x=ixεy=iyεz=izε,ix,iy,iz[Nε/2,Nε/2]Nε=lg(D0/ε)/lg20<ε<<D0
[A1]

then the voxel V0 can be programmably manipulated via its elements indexed by (ix,iy,iz). The Nε×Nε×Nε array for V0 is its support grid. A support grid element takes on a tiny bin in size of ε×ε×ε, which is denoted by Vε. Being the smallest grid element, any property (such as proton spin precession, BOLD field, BOLD signal, etc) associated with a Vε is indivisible and therefore must be manipulated as a whole or by a single entry. In particular, the indivisibility of spin processions of protons in Vε is described as a “spin packet” model in this paper. In principle, the validity of the “spin packet” model demands an as-small-as-possible Vε as long as the support grid array can be managed by a computer program.

Let Vj denote a j-level gridel with the label of (1k,2k,…,jk) in the octree, it occupies a volume space that contains as many as Nj 3 spin packets. For a Vj centered at (ix,iy,iz), we can express its spatial occupancy by

Vj(ix,iy,iz;k1,k2,,kj)={1,ixix<Nj/2,iyiy<Nj/2,iziz<Nj/20,otherwiseNj=round(djε)=round(D0ε2f)=Nε2fcard{Vj(ix,iy,iz;k1,k2,,kj)k1,k2,,kjK8}=8j0<j<<Nε
[A2]

where (ix,iy,iz) denotes the running index for traversing elements in the support grid, and round(.) the rounding operation of taking integer, card{.} counting the number in a set. We see from Eq. [A2] that there are as many as 8j gridels at level j during octadic decomposition.

We can consider the gridel Vj in Eq. [A2] as spatial masks, or mask bins, which follow the spatial containment and disjoint partition relationships, as expressed by

VjVj1V0Vj1(ix,iy,iz;k1,k2,,kj1)=kjK8Vj(ix,iy,iz;k1,k2,,kj1,kj)card{Vj}=8card{Vj1}==8j
[A3]

The cardinality relationship indicates the tree growth by 8j. A J-level octadic partition on voxel V0 produces an octree, as notated by

octreeJ(root=V0)={V0;V1(k1);V2(k1,k2);VJ(k1,k2,,kJ)kjK8,j=1,2,,J}
[A4]

The octree is totally determined by a root and the ultimate resolution level J. It is noted that octreeJ(V0) provides a collection of masks for us to systematically partition a voxel region. For example, we use the multigrid masks to calculate the BOLD fields and signals over multilevel gridels, thereby constructing octrees for multiresolution BOLD fields and signals, as reported in the text.

Appendix B

Multiresolution BOLD signal tree: octreeJ(S0)

Multigrid decomposition upon voxel V0 produces a spatial mask tree, octreeJ(V0) in Eq. [A4]. For each gridel node, we calculate the BOLD signal. Accordingly, we construct an octree to systematically represent the multiresolution analysis of BOLD signal, that is

octreeJ(S0(t))={S0(t),S1(t;k1),S2(t;k1,k2),,SJ(t;k1,k2,,kJ)kjK8,j=1,2,,J}
[B1]

In comparison with octreeJ(V0), which reveals the spatial containment relationships for multigrid voxel partition, octreeJ(S0(t)) depicts the multiresolution spatial decomposition behavior of BOLD signal. Specifically, octreeJ(S0(t)) offers the following properties.

  1. Each node signal is spatially decomposed into 8 subnode signals (indexed by K8), which are governed by an additive identity
    Sj1(t;k1,k2,,kj1)kjK8Sj(t;k1,k2,,kj1,kj),j[1,Nε].
    [B2]
    where Nε represents the maximum resolution level that is determined by the support grid resolution ε (see Eq. [A1]).
  2. The summation of BOLD signals over all nodes at any level j is invariant, that is
    S0(t)k1K8k2k8kjK8Sj(t;k1,k2,,kj),j[1,Nε].
    [B3]
    This property is called global rebinning invariance because all spin packets (in a number of Nj 3) at level j are conceptually rebinned into gridels, but not subject to any nonlinear operation, thus remaining intact.
  3. Any node is invariant to its multilevel decomposition (to any level), as expressed by
    Sj(t;k1,k2,,kj)kj+1K8kj+2K8kj+jK8Sj(t;k1,k2,,k1,k2,,kj,kj+1,,kj+j),j>j
    [B4]
    This property is called local rebinning invariance, as it is applies to a gridel at level j>0.
  4. The amplitudes and phases of gridel BOLD signals at level j follow the inequality relationships
    k1=K8k2=K8kj=K8Aj(t;k1,k2,,kj)=k1=K8k2=K8kj=K8sj(t;k1,k2,,kjA0(t)=S0(t),j[1,Nε]kj+1=K8kj+j=K8Aj(t;k1,k2,,kj,kj+1,,kj+j)=kj+1=K8kj+j=K8sj(t;k1,k2,,kj)Aj(t)=kj+1=K8kj+j=K8Sj(t;k1,k2,,kj),j>j
    [B5]
    and
    k1=K8k2=K8kj=K8Φj(t;k1,k2,,kj)Φ0(t),j[1,N]kj+1=K8kj+j=K8Φj(t;k1,k2,,kj,kj+1,,kj+j)Φj(t;k1,k2,,kj),j>j
    [B6]
    where the equality relationship holds for perfect uniform magnetic field distribution over the gridel in question, implying no contrast.

The properties 1-3 are useful for program verification during algorithm implementation, and the property 4 is responsible for spatial contrast regeneration of magnitude and phase due to multiresolution BOLD signal decomposition.

References

1. Hyde JS, Biswal BB, Jesmanowicz A. High-resolution fMRI using multislice partial k-space GR-EPI with cubic voxels. Magn Reson Med. 2001;46(1):114–25. [PubMed]
2. Kirwan CB, Jones CK, Miller MI, Stark CE. High-resolution fMRI investigation of the medial temporal lobe. Hum Brain Mapp. 2007;28(10):959–66. [PMC free article] [PubMed]
3. Kriegeskorte N, Bandettini P. Analyzing for information, not activation, to exploit high-resolution fMRI. Neuroimage. 2007;38(4):649–62. [PMC free article] [PubMed]
4. Soltysik DA, Hyde JS. High spatial resolution increases the specificity of block-design BOLD fMRI studies of overt vowel production. Neuroimage. 2008;41(2):389–97. [PMC free article] [PubMed]
5. Yoo SS, Guttmann CR, Panych LP. Multiresolution data acquisition and detection in functional MRI. Neuroimage. 2001;14(6):1476–85. [PubMed]
6. Lee J, Shahram M, Schwartzman A, Pauly JM. Complex data analysis in high-resolution SSFP fMRI. Magn Reson Med. 2007;57(5):905–17. [PubMed]
7. Boecker H, Sefat D, Kleinschmidt A. High-resolution functional magnetic resonance imaging of cortical activation during tactile exploration. Human brain mapping. 1995;3(3):161–256.
8. Introduction to functional magnetic resonance imaging: principle and techniques. Buxton RB: Cambridge University Press; 2002.
9. Rowe DB. Modeling both the magnitude and phase of complex-valued fMRI data. Neuroimage. 2005;25(4):1310–24. [PubMed]
10. Rowe DB, Logan BR. Complex fMRI analysis with unrestricted phase is equivalent to a magnitude-only model. Neuroimage. 2005;24(2):603–6. [PubMed]
11. Rowe DB, Meller CP, Hoffmann RG. Characterizing phase-only fMRI data with an angular regression model. J Neurosci Methods. 2007;161(2):331–41. [PMC free article] [PubMed]
12. Haacke EM, Brown R, Thompson M, Venkatesan R. Magnetic resonance imaging physical principles and sequence design. New York: John Wiley & Sons, Inc; 1999. pp. 753–754.
13. Bracewell RN. The Fourier Transform and Its Applications. 3. Boston: McGraw-Hill; 2000.
14. Arja SK, Feng Z, Chen Z, Caprihan A, Kiehl KA, Adali T, Calhoun VD. Changes in fMRI magnitude data and phase data observed in block-design and event-related tasks. Neuroimage 2009 [PMC free article] [PubMed]
15. Feng Z, Caprihan A, Blagoev KB, Calhoun VD. Biophysical modeling of phase changes in BOLD fMRI. Neuroimage. 2009;47(2):540–8. [PubMed]
16. Nan FY, Nowak RD. Generalized likelihood ratio detection for fMRI using complex data. IEEE Trans Med Imaging. 1999;18(4):320–9. [PubMed]
17. Rowe DB. Magnitude and phase signal detection in complex-valued fMRI data. Magn Reson Med. 2009;62(5):1356–7. author reply 1358-60. [PMC free article] [PubMed]
18. Rowe DB, Logan BR. A complex way to compute fMRI activation. Neuroimage. 2004;23(3):1078–92. [PubMed]
19. Xu Y, Wu G, Rowe DB, Ma Y, Zhang R, Xu G, Li SJ. COmplex-Model-Based Estimation of thermal noise for fMRI data in the presence of artifacts. Magn Reson Imaging. 2007;25(7):1079–88. [PMC free article] [PubMed]
20. Zhao F, Jin T, Wang P, Hu X, Kim SG. Sources of phase changes in BOLD and CBV-weighted fMRI. Magn Reson Med. 2007;57(3):520–7. [PubMed]
21. Menon RS. Simulation of BOLD phase and magnitude changes in a voxel. Proc Intl Soc Mag Reson Med. 2003;11:1719.
22. Pathak AP, Ward BD, Schmainda KM. A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: the finite perturber method. Neuroimage. 2008;40(3):1130–43. [PMC free article] [PubMed]
23. Martindale J, Kennerley AJ, Johnston D, Zheng Y, Mayhew JE. Theory and generalization of Monte Carlo models of the BOLD signal source. Magn Reson Med. 2008;59(3):607–18. [PubMed]
24. Chen Z, Molloi S. Automatic 3D vascular tree construction in CT angiography. Comput Med Imaging Graph. 2003;27(6):469–79. [PubMed]
25. Chen Z, Ning R. Supergridded cone-beam reconstruction and its application to point-spread function calculation. Appl Opt. 2005;44(22):4615–24. [PubMed]
26. Chen Z, Ning R. Breast volume denoising and noise characterization by 3D wavelet transform. Comput Med Imaging Graph. 2004;28(5):235–46. [PubMed]
27. Salomir R, de Senneville BD, Moonen CTW. A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts Magn Reson B. 2003;19:26–34.
28. Koch KM, Papademetris X, Rothman DL, de Graaf RA. Rapid calculations of susceptibility-induced magnetostatic field perturbations for in vivo magnetic resonance. Phys Med Biol. 2006;51(24):6381–402. [PubMed]
29. Marques JP, Bowtell R. Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility Concepts Magn. Reson;B. 2005;25:65–78.
30. Kiselev VG, Posse S. Analytical model of susceptibility-induced MR signal dephasing: effect of diffusion in a microvascular network. Magn Reson Med. 1999;41(3):499–509. [PubMed]
31. Hendricks WAR, K W. The sampling distribution of the coefficient of variation. Ann Math Statist. 1936;7(3):4.
32. Rogers T. Amazing applications of probability and statistics. 2010. [March 25, 2010]. http://www.intuitor.com/statistics/CentralLim.html.
33. Boxerman JL, Bandettini PA, Kwong KK, Baker JR, Davis TL, Rosen BR, Weisskoff RM. The intravascular contribution to fMRI signal change: Monte Carlo modeling and diffusion-weighted studies in vivo. Magn Reson Med. 1995;34(1):4–10. [PubMed]
34. Ogawa S, Tank DW, Menon R, Ellermann JM, Kim SG, Merkle H, Ugurbil K. Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc Natl Acad Sci U S A. 1992;89(13):5951–5. [PubMed]
35. Feng Z, Caprihan A, Blagoev KB, C CV. Biophysical Modeling of Phase Changes in BOLD fMRI. NeuroImage. 2009;47(2):540–548. [PubMed]