Absolute Oxygenation Metabolism Measurements Using Magnetic Resonance Imaging

R2 Approaches

The oxygenation dependence of R2 originates from the irreversible dephasing caused by water exchanging and diffusion through the magnetic field gradient induced by deoxyhemoglobins. By using a Carr-Purcell-Meiboom-Gill (CPMG) sequence, the Luz-Meiboom model describes the relationship between R2 and oxygenation levels through the exchange between two sites that have different frequencies as [44]:

where R₂₀ is the R₂ for fully oxygenated blood, τ_ex is the average time required for a proton to move between the two sites at different frequencies; ω₀ is the Larmor frequency, α is a coefficient related to the geometry of the erythrocyte and susceptibility of deoxyhemoglobin; P_A is the fraction of protons resident at one of the sites under exchange; τ₁₈₀ is the interval between the two adjacent π pulses in the sequence; and % HbO2 is the fractional oxygen saturation. The parameters in Eq. (1) depend on factors such as pulse sequence, field strength, and the geometry of the erythrocyte. Some of these parameters, including α, P_A and τ_ex, are difficult to determine. Wright et al [45] proposed an ex vivo calibration method to compute oxygen saturation from this equation without concerning all parameters and a single parameter K, which is a function of τ₁₈₀ and ω₀, is used to absorb all parameters except (1-% HbO2) as

With ex vivo blood samples and in a range of physiologically relevant oxygen saturation, a calibration curve was obtained between R2 of the blood and its oxygen saturation based on Eq. (2) for a specific sequence and field strength. With the assumption that K is identical between the ex vivo and in vivo conditions, this experimentally determined calibration curve was employed for in vivo experiments to obtain oxygen saturation in the aorta, superior vena cava and pulmonary trunk. Results were promising for these vessels and particularly so for venous blood. Moreover, with a similar approach, Flotz et al. [46] demonstrated that a positive correlation existed between MR estimated oxygen saturation at the descending aorta and the readings from blood gas analysis of the blood samples collected from an intra-arterial catheter placed within the descending aorta. Fig (1) demonstrated an overview of the methods. T2 in coronary sinus was computed from a series of spin-echo images with different echo times. An ex vivo calibration curve was then utilized to convert the blood T2 to blood oxygenation.

By assuming that water exchange between erythrocyte and plasma is rapid and diffusion effect is negligible, Van Zijl et al. incorporated previous works [31,47-50] into a more comprehensive theoretical model and derived an expression of blood R2 as a sum of R2 from plasma, erythrocyte, and a relaxation rate due to water exchange [51]:

where R_2,exch is the relaxation rate due to water exchange at two sites with different frequencies [44], i represents arterial, capillary, or venous blood compartment, respectively. With the assumption that water exchange between blood vessels and brain tissue is slow, the signal intensity in the brain parenchyma can be expressed as [51]:

where ${}_{\mathrm{tissue}}+{}_i{}_{\mathrm{blood},i}=1$ R_2,app is the apparent rate measured through fitting the experimental data to a single exponential function. χ_tissue and ${}_i{}_{\mathrm{blood},i}$ are the volume fractions of tissue and blood, respectively. Similar to the Luz-Meiboom equation (Eq 1), multiple factors including geometry, blood volume, pulse sequence, Hct, magnetic field strength, and venous oxygenation can affect R2. After simplyfing this method by eliminating the signal from tissue, Golay et al. and Oja et al. have measured T2 within voxels only containing pure venous blood. Through an ex vivo calibration, they estimated venous blood oxygenation in large cerebral veins during resting and activation states [52, 53].

Along this direction, Lu and Ge developed a method called T2-Relaxation-Under-Spin-Tagging (TRUST) to overcome the limitations of previous approaches through separating blood signal from tissue regardless of partial volume effects [54]. An overview of this method is given in Fig. (2). Spin labeling principle was utilized to preserve venous blood signal while suppressing tissue signal by subtracting ‘label’ from ‘control’ scans. T2 of the isolated blood signal was measured by employing several T2 weighting scans using flow-insensitive T2-preparation pulses [54]. Similar to previous approaches [52, 53], an ex vivo calibration was used to convert T2 into venous blood oxygenation (Fig. 2c). More recently, Xu et al. combined TRUST method with a phase contrast MRI to measure both venous blood oxygen saturation in the superior sagittal sinus and internal jugular veins and blood flow velocity in major arteries, such as the internal carotid artery and vertebral arteries for a calculation of whole brain global CMRO₂. This technique was also utilized to study normal brain aging [55] and the effects of CO₂ on brain oxygen metabolism [56]. Lu et al found that venous oxygen saturation in superior sagittal sinus decreased and the global CBF reduced, resulting in an increased global CMRO₂ as a function of age [55]. Using 5% CO₂ enhanced air, Xu et al demonstrated that venous oxygen saturation and CBF increased, while the global CMRO₂ reduced by 13.4%±2.3% [56]. The TRUST method can single out venous blood without the partial volume effects from tissue. Therefore, this method can be utilized within voxels fully or partially occupied by venous blood. However, it is clear from Eq. (4) that MR signal changes caused by the alteration of R2 rely largely on the available blood volume. Since brain parenchyma has cerebral blood volume (CBV) ranging between 2 to 5% under normal physiological conditions [57-62], it is very challenging to measure brain tissue oxygenation [45, 46, 52-54] because of the low sensitivity of R2. To overcome this hurdle, Bolar et al proposed a voxel based method using venous-targeted velocity selective (VS) spin-labeling method [63]. A cutoff velocity Vc was used to only preserve signal from, presumably, post-capillary venular (PCV) blood accelerated from below Vc to above Vc. Signal from fast flowing arteries (>Vc) was suppressed through velocity selection saturation, while static, CSF and non-venular blood were removed through a subtraction between label and control images. T₂ from the isolated PCV blood were obtained and converted to venous oxygen saturation (Yv) using an ex vivo calibration curve. An estimate of CMRO₂ was then obtained in gray matter using OEF computed from Yv and CBF measured from an ASL scan [63]. This approach, dubbed as QUantitative Imaging of eXtraction of Oxygen and TIssue Consumption (QUIXOTIC) MRI, has been successfully employed to detect venous oxygen saturation increase in response to a visual stimulus [64] and to hypercapnia [65]. To further improve the SNR of QUIXOTIC method, Guo and Wong proposed a strategy to image the PCV blood using VS Excitation with Arterial Nulling (VSEAN) [66]. A slab selective inversion pulse is applied below the imaging plane with an inversion time to null arterial blood. To remove signal from static tissue, two VS modules are utilized to generate a quadratic sine modulation (sin²(πVc)) of longitudinal magnetization. The static and slow moving spins with near zero velocities are suppressed. Compared to QUIXOTIC methods, major advantages of VSEAN method include a higher SNR due to more relaxed venous blood and an improved time efficiency because of no need of a subtraction to remove static tissue signal. Both QUIXOTIC and VSEAN can measure venous oxygen saturation within gray matter but not in white matter due to low signal sensitivity. The cutoff blood velocity assumption may need special attention during altered physiological and pathological conditions.

The advantage of these R2 based approaches lies at their insensitivity to the background magnetic field inhomogeneity. Due to the complex signal dependency on various factors as listed in Eq.(1) – (4), it is very difficult to experimentally determine these factors in vivo. An ex vivo calibration needs to be performed by matching the pulse sequence, imaging parameters, temperature, pH, and Hct in the blood samples with the in vivo scan for accurate calibration. This calibration may be problematic under pathological conditions since many of the physiological parameters may deviate dramatically from their normal values.

R2*/R2' Based Methods

In the presence of paramagnetic materials, spins precess at frequencies different from the Larmor frequency, which is referred to as a frequency shift. Usually, microscopic magnetic susceptibility sources are confined to a specific region or compartment, but their susceptibility effects can extend far beyond the compartment boundary[67]. One such example is the cylindrical body utilized to model blood vessels. Consider a cylinder that has a uniform and constant susceptibility Δχ with respect to background, if the radius, a, is much smaller than the length, the frequency shift outside and inside of the infinitely long cylinder are given in cgs units as [68]:

where θ is the angle between the main magnetic field B₀ and the long axis of a cylinder, a is the cylinder radius, ρ is the distance between the point of interest and the center of the cylinder cross-section in the plane normal to the cylinder, and Ф is the angle between this vector and the component of B₀ in the plane. The signal alteration is a function of susceptibility Δχ as well as the angle θ in Eqs. (5) and (6). In the case of venous blood vessel, $={}_0\mathrm{Hct}\left(1Y\right)$, where Δχ₀ is the susceptibility difference between the fully oxygenated and fully deoxygenated blood and measured to be 0.18 ppm in cgs units [69]; Y is the venous blood oxygen saturation. Hematocrit (Hct) is the volume percentage of RBCs in blood, which ranges between 0.4-0.45 in normal subjects depending on the sex of the subjects [70]. As shown in Eqs. (5) and (6), deoxyhemoglobin behaves as an endogenous contrast agent that cause mesocopic magnetic field variations both within and beyond the venous blood vessels.

Numerical and analytical modeling have been proposed to relate MR signal to tissue oxygenation [49, 71-76]. Two mechanisms affect MR signal formation in the presence of local magnetic field inhomogeneities: static dephasing induced by the phase dispersion of spins experiencing different local magnetic fields and signal alteration arising from the water diffusion through the magnetic field gradient induced by susceptibility perturbers. The static dephasing induced time invariant signal loss can be reversed with a spin echo. In contrast, the incoherent phase dispersion accumulated through water diffusion across the varied magnetic field cannot be recovered. These two effects depend on several factors, such as the absolute value of susceptibility differences, pulse sequences and imaging parameters, size and geometric distribution of the susceptibility perturbers, and possibly the orientation of the sources if the susceptibility pertubers are not symmetric in space. Depending on the relative rates of static dephasing and the water diffusion induced signal loss, changes in MR signal can be separated into three regimes (motion averaged, intermediate and static dephasing) [73]. In motion averaged regime, the diffusion rate is much faster than the static dephasing rate. When the diffusion is slow, the system is said to be in a static dephasing regime. When the diffusion rate is comparable to the static dephasing rate, the system is in an intermediate regime.

A numerical model to simulate the static dephasing regime was proposed by Majumdar in 1991 [71]. Majumdar found that the signal decay in gradient echo images was directly related to the absolute difference of susceptibility and the spatial distribution of the susceptibility sources. Furthermore, he found that at a very short echo time, the gradient echo signal deviated from monoexponential decay in both computer simulation and ex vivo experimental data. Based on Monte Carlo simulations to evaluate both the static dephasing and diffusion effects, Muller et al. [49] and Boxerman et al. [74] demonstrated that transverse relaxation rates R2 and R2* had dependence on vessel size. In gradient echo image, R2* increased along with the increase of perturber size for very small vessels, and it remained unchanged when the perturbers were greater than a certain size. On the other hand, the maximal R2 was obtained in vessels with an intermediate size. The different vessel size dependence of R2* and R2 might reveal the geometric dependence on water proton diffusion effects as well as the interplay between the static dephasing and diffusion effects. These Monte Carlo simulations demonstrated that pulse sequence, echo time, vessel size and concentration of contrast agent had significant effects on relaxation rate, while vascular permeability to water and flow of erythrocytes had little effect on it [74].

Yablonskiy and Haacke proposed an analytical model to characterize the MR signal alteration in the static dephasing regime [72]. Three random geometric distributions of magnetized spheres, parallel magnetized cylinders, and cylinders with random orientations were explored. In brain parenchyma, blood vessels form an interconnecting network of long cylinders (compared to their radius) with random orientations. Therefore, of particular interest is the randomly oriented cylinder model, in which the intravascular signal contribution is not included by assuming a small blood volume fraction and only signal changes from the extravascular compartment are considered. The measured MR signal can be characterized as: (the effects of spin-spin relaxation rate T2 are ignored in the following equations).

where ρ is the spin density; λ is the volume fraction containing deoxyhemoglobin, which is cerebral venous blood volume (vCBV); and $f_c\left(t\right)$ can be described as

where J0(x) is the zeroth order Bessel function; and δω is the characteristic frequency shift. If one assumes that the arterial blood is fully oxygenated then λ represents the venous blood volume and δω is defined as

where Hct is the fractional hematocrit; B0 is the main magnetic field strength; Δχ₀ is the susceptibility difference between fully oxygenated and fully deoxygenated blood which has been measured to be 0.18 ppm per unit Hct in cgs units [69]; and Y is the oxygen extraction fraction. Two asymptotic forms, namely short time scale ($\left(t\right)1.5$ ) and long time scale ($\left(t\right)>1.5$) are given to approximate the signal equation as:

and

where the subscript s and l respectively denote the short and long time scales; tc is the critical time defined as $\mathrm{tc}=1/$. The commonly used transverse relaxation rate R2* in gradient echo images can be approximated as R2*=R2+ R2', with R2 and R2' represent irreversible and reversible signal decay rate, respectively. The reversible signal decay rate of local susceptibilities (R2') can be written as

R2' can be estimated by fitting a linear curve of the logarithm of S_l(t) with t > 1.5/δω. λ can be calculated as $=\ln \left(S_1\left(t=0\right)\right)\ln \left(S_8\left(t=0\right)\right)$, where $S_1\left(t=0\right)$ can be obtained through an extrapolation of Eq. (11) following the estimation of R2'. After both R2' and λ are estimated, measurement of Y can be obtained from Eq. 9.

Similar to the computer simulation results shown by Majumdar in 1991 [71], this analytical signal model [72] shows that the MR signal is a Gaussian function at short echo time and it becomes monoexponential at longer echo time. This signal model established a direct relationship between the susceptibility variations and the MRI measured signal. The close form analytical signal description, which might be approximated by two asymptotic equations, could be utilized to estimate the volume fraction of the susceptibility sources as well as the absolute susceptibility. Using this signal model, Yablonskiy measured volume fraction of susceptibility sources in phantom [77]. Building upon the Yablonskiy and Haacke’s work, Kiselev and Posse [76] developed a more complete analytical signal model by taking into account the signal from the intravascular compartment as well as the effects of slow diffusion. Within a slow diffusion regime, this Kiselev and Posse model agree with a Monte Carlo simulation [74] and previous models [72, 73] concerning diffusion effects. However, it lacks a close form solution that allows for a direct experimental measurement. More recently, based on the previous analytical model as in [72], He and Yablonskiy developed a quantitative BOLD signal model (qBOLD) by including multi-compartmental structures of brain tissue and considering different T2 decay rates in gray and white matters [26]. Moreover, frequency shift between brain tissue and cerebrospinal/interstitial fluid was also considered.

We have utilized the analytical model proposed by Yablonskiy and Haacke to estimate cerebral blood oxygen saturation in normal human subjects under both normal and hypercapnic conditions, induced by inhalation of room air or carbogen (3% CO₂ mixed with 97% O₂), respectively. In consistent with physiological expectation, an increase of vCBV (corresponds to λ in Eq. 7) and a decrease of both R2′ and OEF are obtained under the hypercapnic condition (Fig. 3, lower row) when compared to that under the normocapnic condition (Fig 3, upper row). Across all subjects, OEF decreased from 0.39±0.04 to 0.28±0.02, while vCBV increased from 3.76%±0.78% to 4.4%±0.67% from normocapnic to hypercapnic conditions, in agreement with the previously reported PET results [78-80]. Moreover, a direct comparison between measurements obtained using the MR based methods and the gold standard, blood gas analysis was performed to validate the MR approach in rats over a wide range of physiological relevant blood oxygen saturation [81]. Moderate hypoxia, severe hypoxia, hypercapnia, and normal control states were induced by giving rats two levels of nitrogen/air mixture, carbogen and air to alter brain oxygen saturation fraction from 0.1-0.8 (the full possible range of tissue oxygen saturation is 0-1). As shown in Fig. (4), MR measured oxygen saturation O₂Sat_MR maps showed moderate and severe reduction in cerebral oxygen saturation during moderate hypoxia (1:2 mixture of nitrogen and air), severe hypoxia (1:2 mixture of nitrogen and air) and an increase of cerebral oxygen saturation during hypercapnia (3% CO₂, 97% O₂). Fig. (5) demonstrated that the MR measured cerebral oxygen saturation is in excellent agreement with the gold standard blood gas analysis for a wide range of cerebral oxygenation. Furthermore, utility of this MR based methods have been explored in both human stroke and ischemic rat model [81, 82]. In addition, He and Yablonskiy have utilized the qBOLD model to estimate cerebral oxygen saturation in human and further validated this method by comparing MR measurements with venous blood gas oximetry in rats that under different anesthetic conditions [83]. As shown in Fig. (6), a high level of oxygen saturation was obtained in rats under isoflurane anesthesia (77%) than that under alpha-chloralose anesthesia (67%).

Unlike R2 based methods, R2′ can be directly linked to brain oxygenation without a need of ex vivo signal calibration. In addition, R2′ based methods are more sensitive to the signal alteration induced by paramagnetic deoxyhemoglobin, which makes the quantification of oxygenation within tissue possible. However, R2′ based methods are more subjective to background field inhomogeneity arising from air tissue interfaces, such as in the frontal orbital or lateral temporal areas. We have shown that R2′ and λ will be overestimated if Eq. 7-9 are utilized [84]. In general, background field inhomogeneity is macroscopic in scale (on the order of a few mm or cm) and has a either linear or low-order spatial dependence [85]. In contrast, the scale of magnetic field variation induced by physiological parameters, such as deoxyhemoglobin, is on the order of a few microns and usually smaller than the voxel size and varies randomly [77]. Several methods have been proposed in an attempt to recover the signal loss induced by the macroscopic background field inhomogeneity or to separate the effects of these two field variations. These approaches are to (1) increase spatial resolution [47, 86]; (2) use a tailored RF pulse to compensate for the signal loss [87]; (3) apply z-shimming gradients to recover signal loss [85, 88-91]; (4) employ an image post-processing method to separate the effects of physiological related susceptibility effects from macroscopic background field [77, 92]; and (5) acquire a background field map to correct for its effect [79]. The dimension along z direction is usually larger than the other two directions (x and y) in 2D MRI images, resulting in a larger degree of intravoxel phase dispersion along this direction. The z-shimming method proposed by Yang et al [85] utilized an additional gradient table along z direction to compensate background field variations. Usually, the maximal z-shimming gradient is determined using the magnitude of background field inhomogenieties and TE. If the z shimming gradient table is chosen correctly, one of the z-shimming gradients will move the k space data back to Kz=0. After a FFT along the z direction, the signal loss can be compensated. Fig. (7) demonstrates that z shimming method can effectively reduce background susceptibility artifacts, resulting in a reduced R2* after z shimming correction. To recover signal loss induced by a range of background susceptibility, multiple steps of z-shimming gradient is needed. For example, if an eight step z-shimming table is utilized, an eight times longer data acquisition time is needed without gaining any SNR. Therefore, though z-shimming method is effective in reducing background susceptibility artifacts, it is not time efficient. On the other hand, background magnetic field map can be obtained through a high resolution gradient echo sequence as shown in Fig. (8) (left) [79]. This field map can then be utilized to correct for the overestimated vCBV (Fig. 8). This method assumes that the background field inhomogeneity is linear. If the magnetic field inhomogeneity is large, this method cannot recover the signal loss accurately.

Several assumptions are made for the analytical signal models discussed above [72]. First, numerous blood vessels with random orientation are assumed so that a statistical approach can be utilized to derive the signal equation. This assumption holds true for brain regions consisting of small vessels, but no longer valid for voxels containing a dominant large vessel, such as a pial or a draining vein. Second, only extravascular signal alterations induced by deoxyhemoglobin is considered. Since normal cerebral blood volume ranges between 2-5% in the brain parenchyma [57-62], signal contributions originating from the intravascular space, which is weighted by CBV, are likely to be small. This assumption may not be valid under certain pathophysiological conditions, such as high grade tumors [93] that usually has high rCBV. Appropriate modifications to include intravascular signal contributions will be necessary in this case. Pulse sequences that can minimize the intravascular signal contributions by dephasing flow spin of blood [80] may alleviate this problem. Finally, only static dephasing mechanism is considered in the signal model, while signal reduction caused by incoherent phase shift from diffusion is ignored. It has been suggested that the diffusion effects depend on pulse sequences, TE, and vessel size [74]. Similar to the approach in [74], we performed a Monte Carlo (MC) simulation to evaluate to what extent that diffusion may impact signal and its dependence on vessel size. In this MC simulation, vessels were generated by uniformly sampling a sphere for random vessel orientations, the number of vessels were determined to occupy 2% of total voxel volume. Brownian motion was used to simulate diffusion effects. In total, 1 million spins were simulated. As shown in Fig. (9), our MC simulation clearly demonstrates that diffusion effects causes signal deviation from static dephasing regime for vessels with a diameter of 7 µm, but not for vessels with a diameter of 25 µm. Usually, short echo time and a high magnetic field favor static dephasing assumption [73].

Once tissue oxygen saturation is obtained, OEF can be computed as OEF=O₂Sat_a-O₂Sat_v, where O₂Sat_a is the oxygen saturation in artery which is usually close to 100%. Together with the CBF obtained by using either dynamic susceptibility contrast (DSC) or arterial spin labeling methods [94-97], an oxygen metabolic index can be computed as a product of OEF and CBF [98].