A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO 2 ) with single gas calibrated fMRI

Abstract

One promising approach for mapping CMRO₂ is dual-calibrated functional MRI (dc-fMRI). This method exploits the Fick Principle to combine estimates of CBF from ASL, and OEF derived from BOLD-ASL measurements during arterial O₂ and CO₂ modulations. Multiple gas modulations are required to decouple OEF and deoxyhemoglobin-sensitive blood volume. We propose an alternative single gas calibrated fMRI framework, integrating a model of oxygen transport, that links blood volume and CBF to OEF and creates a mapping between the maximum BOLD signal, CBF and OEF (and CMRO₂). Simulations demonstrated the method’s viability within physiological ranges of mitochondrial oxygen pressure, P_mO₂, and mean capillary transit time. A dc-fMRI experiment, performed on 20 healthy subjects using O₂ and CO₂ challenges, was used to validate the approach. The validation conveyed expected estimates of model parameters (e.g., low P_mO₂), with spatially uniform OEF maps (grey matter, GM, OEF spatial standard deviation ≈ 0.13). GM OEF estimates obtained with hypercapnia calibrated fMRI correlated with dc-fMRI (r = 0.65, p = 2·10⁻³). For 12 subjects, OEF measured with dc-fMRI and the single gas calibration method were correlated with whole-brain OEF derived from phase measures in the superior sagittal sinus (r = 0.58, p = 0.048; r = 0.64, p = 0.025 respectively). Simplified calibrated fMRI using hypercapnia holds promise for clinical application.

Keywords

Calibrated functional magnetic resonance imaging (calibrated fMRI)cerebral metabolic rate of oxygen (CMRO2)hypercapnia hyperoxia oxygen transport modelling

Introduction

Oxidative metabolism provides most of the brain’s energy and is altered in a variety of pathologies such as neurodegenerative and neuroinflammatory diseases, stroke, epilepsy, and migraine.¹ Magnetic resonance imaging (MRI) approaches for mapping baseline (₀) cerebral metabolic rate of oxygen (CMRO_2,0)^2–8 exploit the Fick Principle, that expresses CMRO² as the product of oxygen delivery (the product of oxygen concentration in arterial blood, C_aO_2, and cerebral blood flow, CBF) and oxygen extraction fraction (OEF) measured in either the macrovascular or the microvascular compartment. Macrovascular CBF₀ can be estimated from volume flow rate in large feeding arteries or draining veins using flow encoding sequences,² whereas OEF₀ can be assessed in draining veins using sequences that measure magnetic susceptibility or blood T₂ predominantly affected by the presence of deoxyhemoglobin (dHb).^8,9 Since large vessels feed or drain significant portions of brain, such measures deliver global or regional information at best. Microvascular CBF₀ can be mapped using perfusion-weighted sequences such as Arterial Spin Labelling (ASL). One drawback of ASL is its low contrast in white matter (WM). Moreover, mapping microvascular OEF₀ is challenging because baseline magnetic susceptibility and MR relaxation parameters within a voxel with a small vascular compartment are not uniquely affected by dHb.

Among others,^10,11 dual-calibrated functional MRI (dc-fMRI)^4,12–15 is a promising approach for OEF₀ and CMRO_2,0 mapping. While measuring CBF₀ with ASL, dc-fMRI estimates OEF₀ from the blood oxygen level dependent (BOLD) signal sensitivity to dHb₀. dc-fMRI uses BOLD-ASL recordings, biophysical modelling of BOLD signal¹⁶ and assumed isometabolic hypercapnic and hyperoxic modulations of CBF and C_aO₂ through respiratory stimuli. BOLD sensitivity to dHb₀ is encoded in the maximum BOLD increase, M, corresponding to complete dHb removal. The two respiratory stimuli decouple the contribution to M of OEF₀ and the dHb₀-sensitive cerebral blood volume, CBV_v,0, when Hb concentration in blood [Hb]^4,14 is known. Although dc-fMRI has been applied in exemplar clinical studies,^17–21 its adoption is limited by the low signal to noise ratio (SNR)²² and by the complex gas challenge paradigm required.

Here, we introduce a new calibrated fMRI framework that estimates OEF₀ with only one measurement of M based on one manipulation of brain physiology and a flow-diffusion model of oxygen transport.^13,23,24 The model describes the steady-state oxygen diffusion from the capillaries into the tissue (equal to CMRO₂) as proportional to the product of the mean capillary transit time (mean CTT, MCTT) and the pressure gradient between the capillary bed and the mitochondria (where the proportionality constant is the effective tissue permeability to oxygen, k). Since MCTT can be expressed as the ratio between capillary blood volume (CBV_cap) and CBF, the flow-diffusion model can be incorporated in the formulation of M by substituting CBV_v,0 for an appropriately scaled CBV_cap,0 (with ρ being the scaling factor). This substitution replaces one unknown variable, CBV_v,0, with two unknowns, one proportionality constant, being a function of ρ and k, and the oxygen pressure at the mitochondria (P_mO_2,0). The advantage of the model lies in the tight distributions of the new parameters and on the reduced effect of their variabilities in the estimate of OEF₀, creating a probabilistic mapping of M, C_aO_2,0 and CBF₀ with OEF₀ and CMRO_2,0 as the parameters to be inferred.

This manuscript reports the validation of the novel single gas approach. We term the new approach using a hypercapnic stimulus, hc-fMRI+, and that using a hyperoxic stimulus, ho-fMRI+. The report is divided into four sections. The first section, by exploiting simulations, describes the advantages, the validity, and the robustness to noise of the framework. The second section investigates the new model in-vivo using a dc-fMRI experiment, employing alternating hypercapnic and hyperoxic gas challenges in healthy subjects. We use a global estimate of OEF₀ in the grey matter (GM), obtained with dc-fMRI analysis,²² and we invert the single gas model using only the hypercapnic or the hyperoxic component of the experiment to investigate the distribution of the proportionality constant and P_mO_2,0 across subjects. The third section validates hc-fMRI+ and ho-fMRI+ against dc-fMRI.²² To do so, the two parameters of the model are fixed to the average values obtained from the previous analysis, and the model is inverted to infer OEF₀. Finally, in the fourth section, GM OEF₀ estimates from the different fMRI approaches are compared to whole-brain OEF₀ inferred from a validated MRI sequence performing phase measures in the superior sagittal sinus (SSS) and conventionally termed ‘OxFlow'.^2,25

Methods

Analytical modeling

Here we summarize the analytical model derivation. Please refer to Supplementary Information for a more detailed description.

BOLD model and the Dual-Calibrated fMRI experiment

The rate of signal decay due to dHb, R₂*|_dHb, within a voxel is represented by:^26,27

{R_{2}^{*} |}_{dHb} = A \cdot {CBV}_{v} \cdot {((1 - {S_{v} O}_{2}) \cdot [H b])}^{β}

(1)

where S_vO₂ is venous saturation, [Hb] is the concentration of hemoglobin in blood and CBV_v is the BOLD sensitive blood volume. β (β = 1.3 at 3 T) and A are field strength and vessel geometry dependent constants. For small perturbations of R₂*|_dHb and using the Grubb relation linking fractional changes in CBV_v and CBF, the steady-state BOLD signal can be expressed, within the Davis Model framework, as:^14,28

\frac{Δ BOLD}{{BOLD}_{0}} = T E \cdot A \cdot {CBV}_{v, 0} \cdot {((1 - {S_{v} O}_{2, 0}) \cdot [H b])}^{β} \cdot \{1 - {(\frac{CBF}{{CBF}_{0}})}^{α} {\cdot (\frac{1 - {S_{v} O}_{2}}{1 - {S_{v} O}_{2, 0}})}^{β}\}

(2)

with the maximum BOLD signal M being equal to:

M = T E \cdot A \cdot {CBV}_{v, 0} \cdot {((1 - {S_{v} O}_{2, 0}) \cdot [H b])}^{β}

(3)

The subscript ₀ depicts baseline values, ΔBOLD/BOLD₀ is the relative BOLD signal change_, TE is the sequence echo-time and α is the Grubb exponent (α = 0.38). During an isometabolic manipulation of brain physiology, equation (2) can be expressed as a function of OEF₀ as:

\frac{Δ BOLD}{{BOLD}_{0}}

= T E \cdot A \cdot {CBV}_{v, 0} \cdot {((1 - \frac{C_{a} O_{2, 0}}{φ \cdot [H b]} \cdot (1 - {OEF}_{0})) \cdot [H b])}^{β} \cdot \{1 - {{(\frac{CBF}{{CBF}_{0}})}^{α} \cdot (\frac{1 - \frac{C_{a} O_{2}}{φ \cdot [H b]} \cdot (1 - \frac{{OEF}_{0} \cdot {CBF}_{0} \cdot C_{a} O_{2, 0}}{CBF \cdot C_{a} O_{2}})}{1 - \frac{C_{a 0} O_{2}}{φ \cdot [H b]} \cdot (1 - {OEF}_{0})})}^{β}\}

(4)

with C_aO₂ being the oxygen concentration in arterial blood and φ being the oxygen binding capacity of hemoglobin (φ = 1.34 mL/g). Even when combining together A and CBV_v,0 in equation (4), the equation still has two unknowns making it not possible to solve for OEF₀ through one manipulation of brain physiology. dc-fMRI solves this by performing two independent manipulations: hypercapnia and hyperoxia. However, the approach suffers from low SNR, a problem that has been addressed by regularizing the inversion procedure for OEF₀²⁹ and by using, simulation-trained, machine learning approaches applied to raw recordings.²²

Flow-Diffusion model of oxygen transport

A simple model can be used to describe the steady-state radial oxygen diffusion into the tissue along a straight cylindrical capillary of unit length:²⁴

\frac{d C_{cap} O_{2} (x)}{d x} = - k \cdot T_{cap} \cdot (P_{cap} O_{2} (x) - P_{m} O_{2})

(5)

where C_capO₂ and P_capO₂ are the concentration and the partial pressure of oxygen at a relative position x along the capillary and T_cap is the CTT. k, the effective permeability, combines the effects of the capillary wall and the surrounding brain tissue into a single interface between the plasma and a well-stirred oxygen pool at the mitochondria at end of the diffusion path, at which the pressure of oxygen is equal to

P_{m} O_{2} .

^13,30 CTT in the single straight capillary is then approximated by the MCTT in the capillary bed within the voxel. MCTT is expressed as the ratio between the capillary blood volume (CBV_cap) and CBF. Since P_capO₂ and C_capO₂ quickly equilibrate (less than a few milliseconds), depending upon the nonlinear nature of Hb binding to oxygen described mathematically by the Hill Equation:

{S O}_{2} = \frac{1}{1 + {(\frac{P_{50}}{{P O}_{2}})}^{h}}

(6)

the following can be obtained:

CBF \cdot \frac{d C_{cap} O_{2} (x)}{d x}

= - k \cdot {CBV}_{cap} \cdot (P_{50} \cdot \sqrt[h]{\frac{C_{cap} O_{2} (x)}{{φ \cdot [H b] - C}_{cap} O_{2} (x)}} - P_{m} O_{2})

(7)

where P₅₀ is the oxygen partial pressure when half of Hb is saturated (generally P₅₀ ≈ 26 mmHg; P₅₀ can be inferred from a measure of end-tidal partial pressure of carbon dioxide, P_ETCO₂), and h is the Hill constant (h = 2.8). An approximated closed solution to the differential equation (7) can be made assuming a linear decrease of C_capO₂(x) and an average C_capO₂(x) equal to <C_capO₂(x)>≈φ · [Hb]·(S_aO₂ + S_vO₂)/2 = φ · [Hb] · (1-OEF/2), where S_aO₂ is the arterial oxygen saturation. Integrating equation (7) and equalizing the oxygen loss from the capillary to CMRO₂, the following is obtained:

{CMRO}_{2} = CBF \cdot OEF \cdot {CaO}_{2} = k \cdot {CBV}_{cap} \cdot (P_{50} \cdot \sqrt[h]{\frac{2}{OEF} - 1} - P_{m} O_{2})

(8)

Integration of the Flow-Diffusion model of oxygen transport into the BOLD model for calibrated-fMRI quantification of CMRO₂

CBV_cap is here assumed to be a fraction of CBV_v, i.e., CBV_v = ρ · CBV_cap. Substituting CBV_cap, from equation (8) into equation (4), we obtain:

\frac{Δ BOLD}{{BOLD}_{0}}

= T E \cdot \frac{A \cdot ρ}{K} \cdot \frac{{CBF}_{0} \cdot {OEF}_{0} \cdot {CaO}_{2, 0} \cdot {((1 - \frac{C_{a} O_{2, 0}}{φ [H b]} \cdot (1 - {OEF}_{0})) \cdot [H b])}^{β}}{(P_{50} \cdot \sqrt[h]{\frac{2}{{OEF}_{0}} - 1} - P_{m} O_{2, 0})} \cdot \{1 - {{(\frac{CBF}{{CBF}_{0}})}^{α} \cdot (\frac{1 - \frac{C_{a} O_{2}}{φ [H b]} \cdot (1 - \frac{{OEF}_{0} \cdot {CBF}_{0} \cdot C_{a} O_{2, 0}}{CBF \cdot C_{a} O_{2}})}{1 - \frac{C_{a} O_{2, 0}}{φ [H b]} \cdot (1 - {OEF}_{0})})}^{β}\}

(9)

with the maximum BOLD signal M equal to:

M = T E \cdot \frac{A \cdot ρ}{K} \cdot \frac{{CBF}_{0} \cdot {OEF}_{0} \cdot {CaO}_{2, 0} \cdot {((1 - \frac{C_{a} O_{2, 0}}{φ [H b]} \cdot (1 - {OEF}_{0})) \cdot [H b])}^{β}}{(P_{50} \cdot \sqrt[h]{\frac{2}{{OEF}_{0}} - 1} - P_{m} O_{2, 0})}

(10)

Equations (9) and (10) encode a non-linear mapping of measurable quantities M, C_aO_2,0 and CBF₀ with OEF₀, enabling OEF₀ (and hence CMRO_2,0) to be inferred using a single manipulation of brain physiology. Apart from the constants that can be indirectly inferred (e.g., P₅₀, [Hb]), assumed (e.g., φ, β) or controlled (e.g., TE), the mapping depends on the non-measurable quantities: A, ρ, k and P_mO_2,0. A, having the same origins as β,³¹ can be estimated assuming primarily an extravascular BOLD signal and assuming R₂*|_dHb $=$ R₂′^32,33. With an experimentally determined cortical R₂′ of approximately 3 s⁻¹ at 3 T³⁴, an average [Hb] of 14 g/dL, a S_vO₂ of 0.6, and a mean CBV_v of 2.5%, from equation (1) we expect a value of A ≈ 14 s⁻¹g^−βdL^β at 3 T. In-vivo variation in ρ has not been studied directly; we discuss this in the Supplementary Information. We expect ρ to be in the range 2 to 3, assuming a capillary blood volume between 20% to 40% of total blood volume, when the arterial contribution is assumed to be 20% to 30%.³⁵ Moreover, we expect a value for the oxygen effective permeability k of around 3 μmol/mmHg/ml/min.²² This value is derived from the literature using a different formalism where oxygen diffusion is assumed to happen at the endothelial wall of capillaries.³⁶ In equations (9) and (10), we create a practical grouping of A, ρ and k into one multiplicative parameter A · ρ/k. At a fixed field strength, all the three parameters are related to tissue structure and vessel geometry, which plausibly affects water and oxygen diffusion in the intravascular and extravascular spaces as well as the volumetric relationship between capillaries, venules and veins. We expect a value of A · ρ/k of the order of A · ρ/k ≈ 10 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min). The mitochondrial oxygen partial pressure at rest, P_mO_2,0, must lie between 0 mmHg and the average oxygen tension of the capillary bed. Several in vivo studies suggest that oxygen tension at brain mitochondria is small in the healthy brain,^23,37 and this theory is consistent with functional hyperemia in response to increased brain oxygen demand. However, departure from a negligible oxygen tension is plausible in the diseased brain.

In summary, the non-linear mapping in equation (9) permits estimation of OEF₀ from one manipulation of brain physiology. Uncertainty in the mapping is driven by variability in two non-measurable quantities, a proportionality constant A · ρ/k, that depends on tissue and micro-vessel structure at a fixed field strength, and P_mO_2,0. Importantly, these non-measurable quantities affect the non-linear mapping differently. The advantage of the new framework lies in the low variability of these parameters and their diminished influence on the OEF₀ estimation compared to CBV_v,0.

Simulations

We performed simulations to investigate the ability of hc-fMRI+ and ho-fMRI+ to infer OEF₀. A forward model using equation (9) was implemented to simulate the BOLD signal and was inverted to retrieve OEF₀. In the forward model some variables were fixed (TE, α, β, h ε, φ) (ε is the oxygen plasma solubility, ε = 0.0031 mL/mmHg/dL) while others were simulated based on random sampling from physiologically and physically plausible distributions. When inverting the model, some random variables were unknown and were either fixed a-priori (A · ρ/k, P_mO_2,0), or inferred (OEF₀, CBV_cap,0 and MCTT_,0). Firstly, we ran the full forward and inverse analysis without measurement noise as a function of either the value chosen a-priori for the random variables that were fixed during the inversion or other parameters of interest (P_mO_2,0 and MCTT_,0). Secondly, we evaluated the effect of measurement noise, which was introduced on measures with lower signal to noise ratio (SNR), namely ASL CBF/CBF₀ and ΔBOLD/BOLD₀. 10⁷ simulations per condition were conducted; the non-linear inversions were performed through explicit search of OEF₀ that explained the measures. The explicit search was performed in the full OEF₀ space (between 0 and 1) with a resolution of 0.01. Constant parameters were set to α = 0.38, β = 1.3, h = 2.8, ε = 0.0031 mL/mmHg/dL, φ = 1.34 mL/g, TE = 30 ms, whereas random variables were simulated using either normal (N) or gamma (Γ) distributions; additional physiological constraints were applied (please refer to the Table in Supplementary Information for additional information).

Figure 1 reports the distributions of the main random variables used in the forward model simulations. With respect to the parameters that were not measured for the inversion, A · ρ/k was simulated using a normal distribution with an average value of 10 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) and a coefficient of variation (CoV) of 0.3, whereas P_mO_2,0 was simulated using a gamma distribution, allowing variation between zero and <P_cap,0> to simulate a large variability in P_mO_2,0 that might be present in disease.

Figure 1.

Random variables used to simulate BOLD and ASL signals using a hc-fMRI+ or ho-fMRI+ forward modelling framework. The variables reported in light grey were assumed to be measured for the hc-fMRI+ or ho-fMRI+ inversion model, those reported in medium grey were fixed a-priori in the inversion model and those in dark grey were inferred by the inversion model.

MRI experiment

Twenty healthy volunteers (13 males, mean age 31.9 ± 6.5 years) were recruited at CUBRIC, Cardiff University, Cardiff, UK. The study was in accordance with the Declaration of Helsinki and was approved by the Cardiff University, School of Psychology Ethics Committee and NHS Research Ethics Committee, Wales, UK. Written consent was obtained from each participant. Data were acquired using a Siemens MAGNETOM Prisma (Siemens Healthcare GmbH, Erlangen) 3 T clinical scanner with a 32-channel receiver head coil (Siemens Healthcare GmbH, Erlangen). A 18 minutes dc-fMRI scan was acquired with interleaved periods of hypercapnia, hyperoxia and medical air being delivered according to the protocol previously proposed.^13,29 3 periods of hypercapnic gas challenges and 2 periods of hyperoxic gas challenges were performed. CO₂ and O₂ in the lungs were evaluated from the volunteer's facemask using a gas analyzer (AEI Technologies, Pittsburgh, PA, USA).

Calibrated fMRI data were acquired during the gas challenge scheme using a pCASL acquisition with pre-saturation and background suppression³⁸ and a dual-excitation (DEXI) readout.³⁹ The labelling duration (τ) and the Post Label Delay (PLD) were both set to 1.5 s, GRAPPA acceleration (factor = 3) was used with TE₁ = 10 ms and TE₂ = 30 ms. An effective TR of 4.4 s was used to acquire 15 slices, in-plane resolution 3.4 mm × 3.4 mm and slice thickness 7 mm with a 20% slice gap. A calibration (S₀) image was acquired for ASL quantification with pCASL labelling and background suppression pulses switched off, with TR = 6 s, and TE =10 ms.¹³ A high-resolution whole brain structural image, used for GM identification in the fMRI space, was acquired using a 3D Fast Spoiled Gradient-Recalled-Echo T1-weighted acquisition (resolution = 1 × 1 × 1 mm³, TE = 3.0 ms, TR = 7.8 ms, TI = 450 ms, flip angle = 20°).

For susceptometry-based oximetry, a transverse slice was acquired at approximately 15 mm above the confluence of sinuses (location at which the inferior sagittal, straight, and transverse sinuses join the SSS) using a T2*-weighted spoiled multi-echo gradient-recalled echo (GRE) sequence with: in-plane resolution = 1.6 × 1.6 mm², slice thickness = 5 mm, field of view (FOV) = 208 × 208 mm², bandwidth= 260 Hz/pixel, three echo times (TEs = 3.92, 7.44, and 10.96 ms), bipolar gradient readout, TR = 35 ms, flip angle = 25°, and acquisition time = 1 min and 7 s. This acquisition was performed in the framework of the OxFlow method, which previous studies have described in detail.^40–42 For vessel identification purposes, two-dimensional T2*-weighted time-of-flight (TOF) images were acquired using a spoiled GRE sequence with: in-plane resolution = 0.86 × 0.86 mm², slice thickness = 2 mm, slice gap = 1.34 mm, FOV = 219 × 219 ×234 mm³, in-plane acceleration factor = 2, bandwidth = 220 Hz/pixel, TE = 4.99 ms, TR = 20 ms, flip angle = 60°. Blood samples were drawn via a finger prick before scanning and were analyzed with the HemoCue Hb 301 System (HemoCue, Ängelholm, Sweden) to calculate [Hb].

fMRI data processing

Gas recordings processing

P_ETCO₂ and P_ETO₂ were extracted from CO₂ and O₂ recordings using in-house software in Matlab (Mathworks, Natick, MA). P_ETCO₂ and P_ETO₂ points were interpolated (cubic spline function), resampled to match fMRI, and shifted in time to maximally correlate with fMRI signals. P_ETCO_2,0 and P_ETO_2,0, were evaluated at baseline in the first 110 seconds. P_ETO₂ was assumed equal to PaO₂ for C_aO₂ computation whereas P_ETCO₂ was assumed equal to PaCO_2. P₅₀ was inferred from estimates of resting blood pH based on the Henderson-Hasselbalch Equation, assuming [HCO₃⁻]= 24 mmol/L:⁴³

p H = 6.1 + \log (\frac{{[HCO}_{3}^{-}]}{0.03 \cdot {P_{a} C O}_{2}})

(11)

and calculating P₅₀ according to the linear relation, P₅₀ = 221.87–26.37 · pH¹³.

SaO₂ was calculated from PaO₂ using equation (6) and CaO₂ was inferred using the relation:

{CaO}_{2} = φ \cdot [H b] \cdot {S O}_{2} + ε {\cdot P O}_{2}

(12)

Finally, to highlight hypercapnic and hyperoxic modulations, P_ETCO₂ and P_ETO₂ traces were high-pass filtered with a 4th order Butterworth digital filter and a high-pass frequency of 1/600 Hz.

fMRI processing

Both functional and structural MRIs were processed using FSL⁴⁴ and in-house algorithms implemented in Matlab. fMRI timecourses were motion corrected based on 6 degrees of freedom co-registration using MCFLIRT.⁴⁵

High-resolution structural T1-weighted MRIs were skull-stripped using BET⁴⁶ and probability maps of Cerebrospinal Fluid (CSF), WM and GM, were computed using FAST.⁴⁷ Motion-corrected fMRI timecourses and the skull-stripped T1-weighted MRI, together with tissue probability maps, were coregistered, relying on 12 degrees of freedom affine transformation, to the S₀ image.⁴⁵ ASL control-tag difference perfusion data (ΔS) in S₀ space were obtained through surround subtraction of the fMRI timecourses at TE₁, normalized with respect to S₀ and converted to CBF in quantitative units of ml/100g/min through the pCASL single compartment kinetic model of labelled spins and voxelwise signal normalization:⁴⁸

CBF = \frac{6000 \cdot λ \cdot e^{\frac{PLD}{{T 1}_{b}}}}{η \cdot η_{inv} \cdot {T 1}_{b} \cdot (1 - e^{- \frac{τ}{{T 1}_{b}}})} \cdot (\frac{Δ S}{S_{0}})

(13)

where λ is the water partition coefficient (λ = 0.9 mL/g), T1_b is the T1 relaxation constant of blood, η is the tagging inversion efficiency (η = 0.85), and η_inv is a scaling factor to account for the reduction in tagging efficiency due to background suppression (η_inv = 0.88).⁴⁹ The T1_b was calculated from SaO₂ and PaO₂ measures using the experimental relation presented in:⁵⁰

{T 1}_{b} = \frac{1}{1.527 \cdot 10^{- 4} \cdot {P_{a} O}_{2} + 0.1713 \cdot (1 - {S_{a} O}_{2}) + 0.5848}

(14)

CBF₀ was evaluated in the first 110 seconds. Finally, fractional CBF was high pass filtered with a 4th order Butterworth digital filter with a high-pass frequency of 1/600 Hz. BOLD T2*-weighted time-courses were obtained through surround averaging of the fMRI at TE₂ and they were expressed as relative BOLD changes with respect to the temporal average of the BOLD signal in the first 110 seconds (BOLD₀). BOLD relative changes were high pass filtered with a 4th order Butterworth digital filter with a high-pass frequency of 1/600 Hz.

Both processed CBF and BOLD volumes were masked with a GM mask at 50% probability threshold.

Dual-calibrated fMRI analysis

Firstly, OEF₀ maps were obtained with a dc-fMRI analysis. Because of the method’s known low SNR, explicit inversion methodologies were avoided and a state-of-the-art method to analyze the data relying on a machine learning approach was used. The machine learning algorithm was fed with fMRI timecourses and, through a time-frequency transformation of fMRI signals to extract features of interest, directly mapped OEF₀ and CMRO_2,0 relying on a pre-trained model based on simulated data. Please refer to²² for detailed information.

Single gas calibrated fMRI analysis

Single gas calibrated fMRI analysis was performed on either the hypercapnic (using hc-fMRI+) or the hyperoxic (using ho-fMRI+) modulations within the dc-fMRI experiment. The evaluation of BOLD and ASL changes with physiological manipulations was performed using the general linear model (GLM).⁵¹ P_ETCO₂ and P_ETO₂ were concurrently regressed on BOLD and ASL filtered modulations. The GLM β-weight delivered an estimate of BOLD or CBF modulation per unit of mmHg of P_ETCO₂ and P_ETO₂. The total modulation was then obtained by multiplying the β-weight with the maximum P_ETCO₂ or P_ETO₂ modulation. The SNR of the modulation was estimated by dividing the GLM β-weight by its confidence interval. In-vivo data were used to evaluate the between subjects distribution of the unknown parameters of the extended model, namely A · ρ/k and P_mO_2,0. This analysis was performed by extracting average BOLD and ASL modulations in the GM. These average estimates were used, together with a global estimate of GM OEF₀ obtained with the dc-fMRI analysis, to invert the model and estimate the unknown parameters. The inversion relied on equation (9), which clearly could not be solved for the two unknowns; however, since A · ρ/k and P_mO_2,0 differently affect the non-linear mapping between BOLD and ASL modulations and OEF₀, we were able to get insight into the average value of both parameters. In particular, we inverted the model assessing the proportionality constant A · ρ/k as a function of the a-priori fixed P_mO_2,0. We expected the A · ρ/k distribution to have a smaller CoV when P_mO_2,0 was closer to the correct average value. The non-linear inversion was performed through an explicit search in the range, for A · ρ/k, between 0 and 40 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) with a resolution of 0.2 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) and, for P_mO_2,0, between 0 and 50 mmHg with a resolution of 1 mmHg. When making predictions with the single gas calibrated models, the unknown parameters were fixed both spatially and between subjects to the optimal values derived in the first step. Using these fixed parameters hc-fMRI+ and ho-fMRI+ inversion models were used for voxelwise estimation of OEF₀ and CMRO_2,0 and comparison with the estimates derived from the dc-fMRI analysis were made.

OxFlow data processing

For a subset of twelve subjects, GM estimates of OEF₀ using single or dual calibrated fMRI approaches were compared to whole-brain estimates of OEF₀ from SSS derived using the OxFlow procedure. OxFlow images were processed using Matlab and code developed in-house. OEF₀ measurements were obtained based on the normalized difference in signal phase between the first and third TEs (Δϕ/ΔTE), with acquisitions having equal gradient polarity.⁴¹ The static background field inhomogeneity was removed using a second-order polynomial fitting.⁴² The intravascular phase was measured as the average signal phase in a region of interest centered in the cross-section of SSS relative to the average signal phase in the tissue region surrounding the SSS. The angle (θ) between the SSS and B₀ was evaluated by comparing the slice acquired for OxFlow and the SSS orientation in the slices immediately above and immediately below in the TOF image. Individual measurements of hematocrit (Hct, %) were obtained based on [Hb] assuming a ratio Hct/[Hb]=3 (%dL/g).⁵²

OEF₀ was calculated using the infinite cylinder analytical model:⁴¹

{OEF}_{0} = \frac{\frac{2 Δ ϕ}{Δ T E}}{γ \cdot Δ χ_{d o} \cdot Hct {\cdot B}_{0} \cdot (\cos^{2} θ - \frac{1}{3})}

(15)

where γ is the proton gyromagnetic ratio (γ = 267.52 · 10⁶ rad/s/T), and Δχ_do = 4π · 0.27 · 10⁻⁶ is the magnetic susceptibility difference between fully oxygenated and fully deoxygenated red blood cells.⁵³

Statistical analysis

Pearson’s correlations and t-tests were performed to assess pairwise associations and biases between the different estimates. Null-hypothesis probabilities (p-values) were calculated using the Student's t distribution (using transformation of correlation for association testing). Normality evaluation was performed prior to statistical inference using the Kolmogorov-Smirnov test.

Results

Simulations

Figure 2 reports the outcome in estimating OEF₀ when using hc-fMRI+ and ho-fMRI+ inversion models with fixed a-priori parameters. Figure 2(a) displays the OEF₀ root mean square error (RMSE) obtained for hc-fMRI+ and ho-fMRI+ with A · ρ/k = 10 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) as a function of P_mO_2,0. A minimum RMSE of OEF₀ = 0.039 was obtained for hc-fMRI+ and a minimum RMSE of OEF₀ = 0.051 was obtained for ho-fMRI+, both at P_mO_2,0 = 11 mmHg. Figure 2(b) displays the scatterplots of the simulated OEF₀ vs. the estimated OEF₀ for hc-fMRI+ and ho-fMRI+ when marginalizing the other variables. The scatterplots reported were obtained using a close to optimal P_mO_2,0, P_mO_2,0 = 10 mmHg, and P_mO_2,0 = 0 mmHg. Figure 2(c) reports the OEF₀ RMSE for the two methods evaluated as a function of two physiological parameters of interest in the forward model, namely MCTT₀ and P_mO_2,0, when fixing a-priori the non-measurable parameters analogous to Figure 2(b). Importantly, when adding noise to BOLD and fractional changes in CBF, the analysis highlighted the stability of the approach with respect to measurement SNR, with the OEF₀ RMSE reaching the OEF₀ RMSE related to model parameters uncertainty at BOLD and ASL SNRs around 4 (refer to Supplementary information for additional information).

Figure 2.

(a) RMSE in OEF₀ for hc-fMRI+ and ho-fMRI+ inversion models with A·ρ/k = 10 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) as a function of P_mO_2,0. (b) Scatterplots of the simulated and estimated OEF₀ for hc-fMRI+ and ho-fMRI+ inversion models assuming either P_mO_2,0 = 0 mmHg or P_mO_2,0 = 10 mmHg; (c) RMSE in OEF₀ as a function of the forward model MCTT₀ and P_mO_2,0 for hc-fMRI+ and ho-fMRI+ inversion models assuming either P_mO_2,0 = 0 mmHg or P_mO_2,0 = 10 mmHg. ***p < 10⁻³.

In-Vivo evaluation of gas, CBF and BOLD modulations

Figure 3 reports the processing steps, in an exemplar subject, that were used to derive the hypercapnic and the hyperoxic CBF and BOLD modulations. Figure 3(a) shows O₂ signals acquired through the gas analyzer with the estimated P_ETO₂ traces whereas Figure 3(b) depicts CO₂ signals and P_ETCO₂ traces. Figure 3(c) shows example of the ASL CBF/CBF₀ and the filtered and fitted P_ETCO₂. Figure 3(d) shows the relative BOLD change and the filtered P_ETCO₂ and P_ETO₂ traces fitted onto ΔBOLD/BOLD₀. The GLM β-weight (in units of cerebrovascular reactivity, CVR, or in units of signal per mmHg of P_ETO₂) were multiplied by the maximum gas modulation to obtain the hypercapnic CBF/CBF₀ and the hypercapnic as well as hyperoxic ΔBOLD/BOLD₀ modulations. Additional information on gas and signal modulations are reported in the Supplementary Information. The GLM analysis delivered an SNR (evaluated as the statistical relevance of the β-weight) of SNR_ASL = 6.1 (SD = 5.4), hypercapnic SNR_BOLD = 16 (SD = 12.7) and hyperoxic SNR_BOLD = 8.6 (SD = 7.52).

Figure 3.

Example of: (a) O₂ and estimated P_ETO₂ traces; (b) CO₂ and estimated P_ETCO₂ traces; (c) GM CBF/CBF₀ and fitted P_ETCO₂ trace. The β-weight of the GLM fit, with units of a CVR, CBF/CBF₀/mmHg, was multiplied by the maximum modulation ΔP_ETCO₂ to obtain the hypercapnic CBF/CBF₀. (d) GM average ΔBOLD/BOLD₀ and fitted P_ETCO₂ and P_ETO₂ traces. The β-weights, with units of %BOLD/mmHg of P_ETCO₂ and P_ETO₂, were multiplied by the maximum modulation ΔP_ETO₂ and ΔP_ETO₂ to obtain the hypercapnic and the hyperoxic ΔBOLD/BOLD₀.

In-Vivo estimation of modeling parameters

Figure 4 reports the analysis performed in-vivo to evaluate the modelling parameters. Figure 4(a) reports the subjects’ average value (and standard error, SE) of A · ρ/k as a function of P_mO_2,0. The value is reported for both hc-fMRI+ and ho-fMRI+. In agreement with equation (9), for higher P_mO_2,0 the estimate of A · ρ/k decreased. We obtained, for a P_mO_2,0 = 0, an average value of A · ρ/k = 8.85 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) (SE = 0.58 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min)) for hc-fMRI+ and A · ρ/k = 6.03 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) (SE = 0.41 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min)) for ho-fMRI+. Figure 4(b) reports the CoV of A · ρ/k for hc-fMRI+ and ho-fMRI+ as a function of P_mO_2,0. The smallest CoV was obtained with a P_mO_2,0 ≈ 0 for both hc-fMRI+ (CoV = 0.29) and ho-fMRI+ (CoV = 0.31) with a monotonic CoV increase at increasing P_mO_2,0.

Figure 4.

Results of the analysis evaluating the modelling unknown parameters that used hc-fMRI+ or ho-fMRI+ and the OEF₀ derived from the dc-fMRI analysis. (a) Subjects’ average (and SE) estimate of the scaling parameter A·ρ/k of the model as a function of the P_mO_2,0 assumed. (b) Subjects’ CoV of the scaling parameter A·ρ/k as a function of P_mO_2,0 assumed. (c) Comparison between hc-fMRI+ and ho-fMRI+ estimates of A·ρ/k for each subject, assuming a P_mO_2,0 = 0 mmHg. **p < 0.01.

Figure 4(c) reports the comparison between hc-fMRI+ and ho-fMRI+ estimates of A · ρ/k for each subject, when fixing the P_mO_2,0 at the value of P_mO_2,0 = 0 mmHg . A good correlation was obtained with a r = 0.71, df = 18, p = 4.2 · 10⁻⁴, with a smaller hyperoxic estimate.

In-vivo estimation of oxygen extraction fraction: Calibrated fMRI vs. dual-calibrated fMRI (section 3)

Figure 5(a) reports exemplar OEF₀ and CMRO_2,0 maps obtained with dc-fMRI, hc-fMRI+ and ho-fMRI+.

Figure 5.

(a) Exemplar GM OEF₀ and CMRO_2,0 maps for a participant of the study obtained with dc-fMRI (left colum), hc-fMRI+ (central column) and ho-fMRI+ (right column). (b) Scatterplots and Bland-Altmann plots comparing the average OEF₀ in the GM between the hc-fMRI+ (upper row) and ho-fmri+ (lower row) and dc-fMRI. **p < 0.01; ***p < 10⁻³.

Notably, subjects’ average spatial variabilities (estimated as standard deviation, SD) in the GM OEF₀ of SD = 0.17 (SE = 0.003), SD = 0.13 (SE = 0.002) and SD = 0.15 (SE = 0.002) were obtained for dc-fMRI, hc-fMRI+ and ho-fMRI+, respectively.

Figure 5(b) reports the scatterplots and the Bland-Altmann plots comparing the average OEF₀ in the GM between dc-fMRI and the single calibration approaches. Average global GM OEF₀ (mean±SD) were 0.39 ± 0.04, 0.39 ± 0.03, and 0.40 ± 0.03 for dc-fMRI, hc-fMRI+ and ho-fMRI+, respectively. hc-fMRI+ OEF₀ was significantly correlated with that of dc-fMRI (r = 0.65, df = 18, p = 2 · 10⁻³, OEF₀ RMSE = 0.033) whereas that of ho-fMRI+ was not (r = 0.26, df = 18, p = 0.27, OEF₀ RMSE = 0.044). No significant bias between the different approaches was found, but this was dependent on the proportionality constant calibration using dc-fMRI. A significant correlation, with no bias, was obtained between hc-fMRI+ and ho-fMRI+ (r = 0.50, df = 18, p = 0.02).

In-vivo estimation of oxygen extraction fraction: Calibrated fMRI vs. OxFlow (section 4)

Figure 6 reports the scatterplots and the Bland-Altmann plots comparing the global OEF₀ of the fMRI approaches to the OEF₀ estimated in the SSS using OxFlow in a subset of 12 subjects. Figure 6(a) shows, for one representative subject, the magnitude image and the processed phase image used to estimate OEF₀ in the SSS within OxFlow. Average SSS OEF₀ estimated using OxFlow was 0.31 ± 0.07. Significant associations of the average OEF₀ in the GM using a fMRI approach with whole-brain OEF₀ retrieved using OxFlow were obtained for dc-fMRI (r = 0.58, df = 10, p = 0.048, RMSE = 0.034, Figure 6(b)) and hc-fMRI+ (r = 0.64, df = 10, p = 0.025, Figure 6(c), RMSE = 0.041). No significant association was obtained using ho-fMRI+ (r = 0.36, df = 10, p = 0.24, Figure 6(d), RMSE = 0.066). A systematic bias was obtained with the OxFlow underestimating the OEF₀ with respect to fMRI. For the two fMRI approaches that delivered a significant association with OxFlow, an absolute difference between the dc-fMRI and OxFlow of ΔOEF₀ = 0.077 with a t = 4.57, df = 11, p = 8 · 10⁻⁴, and an absolute difference between hc-fMRI+ and OxFlow OEF₀ of ΔOEF₀ = 0.083 with a t = 5.13, df = 11, p = 3.2 · 10⁻⁴ were obtained.

Figure 6.

Scatterplots and Bland-Altmann plots comparing the OEF₀ of the calibrated fMRI approaches and OxFlow in a subset of subjects. (a) Example of magnitude (arbitrary units) and processed phase images used to estimate OEF₀ in the SSS within the OxFlow method. SSS and the reference region are outlined in blue and yellow, respectively. OxFlow vs. (b) dc-fMRI; (c) hc-fMRI+; (d) ho-fMRI+. *p < 0.05.

Discussion

We introduced a framework for mapping OEF₀ and CMRO_2,0 using single gas calibrated fMRI. The method integrates a flow-diffusion model of oxygen transport²⁴ with the steady-state BOLD signal model.¹⁶ Simulations suggest the approach to be valid over a wide range of brain physiology. The new approach, when applied to hypercapnia, compared well with dc-fMRI and whole-brain OEF₀ assessed in the SSS using OxFlow. Compared to dc-fMRI, the novel method permits a simpler stimulation paradigm based on a single exogenous gas challenge²⁷ or, presumably, on an endogenous challenge such as breath hold,⁵⁴ and makes the approach robust to measurement noise.

Simulations

The simulations relied on a forward model assuming the new framework to be correct. When inverting the model, the unknown random variables were: (i) A · ρ/k, a lumped parameter dependent on field strength, tissue structure and vessel geometry and (ii) the mitochondrial oxygen pressure at rest, P_mO_2,0. Variability in A · ρ/k (CoV = 0.3) was based on in-vivo data (Figure 4), whereas the P_mO_2,0 was simulated in the range 0-<P_capO₂> (Figure 1). When inverting the forward model with these parameters fixed, we obtained low OEF₀ RMSE (around 0.05) for both hc-fMRI+ and ho-fMRI+ when marginalizing all other variables, with slightly better performance for hc-fMRI+ (Figure 2(a) and (b)). This highlighted the unknown parameters’ reduced effect on the mapping between the measurable variables and OEF₀. In fact, when considering OEF₀ RMSE as a function of a wide range of two interesting physiological variables, P_mO_2,0 and MCTT₀, the OEF₀ RMSE was small. Only with very high P_mO_2,0 (>25 mmHg) and long MCTT₀ (>2.5 s) the RMSE increased significantly. Very high P_mO_2,0 and long MCTT₀, associated with low CBF₀, are expected only in diseases that heavily alter oxygen supply, vasculature and mitochondrial function.

The simulations revealed the effect of BOLD and ASL measurement noise. For both BOLD and ASL modulations, the OEF₀ RMSE quickly reached the value caused by uncertainty in physiology at an SNR ≈ 4. This is the SNR of the modulation estimate, not the temporal SNR of the raw signals. For example, when the modulation is estimated within a GLM framework regressing P_ETO₂ and P_ETCO₂ onto BOLD and ASL modulations (Figure 3), the SNR is the GLM β-weight divided by its confidence interval. Average voxel SNRs were between 6 and 16 in vivo for both signals and gas challenges. The robustness to noise of the approach is advantageous compared to dc-fMRI, that often relies on constrained inversion algorithms, trading off accuracy for higher stability.²²

Modeling parameters

Investigation of model parameters suggested an average value of A · ρ/k of the order of 10 s⁻¹g^−βdL^β/(μmol/mmHg/ml/min) when using hc-fMRI+, and an average value of P_mO_2,0 in the healthy population close to 0 for both hc-fMRI+ and ho-fMRI+ (Figure 4) both results agreed with expectations.^36,55 In fact, the estimate of A · ρ/k decreased beyond expectations at increasing P_mO_2,0 and increased its CoV as a function of the assumed P_mO_2,0, with a rapid increase above 20 mmHg. This work indeed suggests a particularly low average PmO2_,0 in the healthy brain. However, it should be stressed that a strong increase in CoV of A · ρ/k was only observed at PmO_2,0 above 20 mmHg. The confidence interval of the estimate still cannot provide a definitive answer on the average PmO_2,0 within the range 0-20 mmHg. There is work suggesting the mitochondrial PmO_2,0 is about ∼12 mmHg^56,57 which goes against the common assumption of PmO_2,0 being near zero in the healthy brain³⁶ and, indeed, this is still an open debate. Nonetheless, the simulations of the study clearly demonstrate that the approach estimating OEF₀ has limited sensitivity to the value of PmO_2,0 if PmO_2,0 and MCTT₀ are not both very high. The good correlation between the hypercapnic and the hyperoxic estimates indicated consistency. However, we obtained a value of A · ρ/k for ho-fMRI+ around 30% smaller than expected. This result might be a cross-talk effect of the hypercapnic on the hyperoxic BOLD modulations in the dc-fMRI experiment, or an overestimation of ΔPaO₂.

Comparison with dc-fMRI and OxFlow

Spatial homogeneity of OEF₀ in healthy subjects is often taken as an indicator of successful OEF₀ mapping. hc-fMRI+ and ho-fMRI+ decreased OEF₀ spatial variability in GM compared to dc-fmri. The lower variability of hc-fMRI+ and ho-fMRI+ suggests a greater robustness, with respect to measurement SNR, compared to the dc-fMRI. Comparison of GM OEF₀ estimates suggests that hc-fMRI+ is a valid alternative to dc-fMRI (Figure 5(b)). In addition, when comparing GM OEF₀ of the different fMRI approaches with global OEF₀ in the SSS through OxFlow, clear associations with the OxFlow OEF₀ were obtained for both dc-fMRI and hc-fMRI+ (Figure 6). We identified a bias between the OxFlow and the fMRI estimates. OxFlow yielded a lower global OEF₀ with values around 0.31 and 0.39 for Oxflow and fMRI, respectively. This difference might be explained by previous work, where Oxflow using analytical modelling gave high estimates of venous saturation.⁵⁸ Systematic differences in the global OEF₀ estimated using heterogeneous MRI or non MRI techniques are established in the literature. Work performed using measures of T2 relaxation in venous blood, delivered global OEF₀ between 0.36 and 0.38.^59,60 These values were close to global OEF₀ estimated using the gold standard ¹⁵O positron emission tomography (¹⁵O PET)^60,61 and are more compatible with the OEF₀ obtained with calibrated fMRI approaches. Nonetheless, once a good correlation between modalities is established, biases can be corrected using data-driven approaches or through better models of the underlying physics and physiology. The low performance of ho-fMRI+ is indeed a negative result of the study. Hyperoxia is generally better tolerated than hypercapnia¹⁹ and it would be more easily applicable in clinical settings. The lower performance of hyperoxia is plausibly related to the noisier estimate of M. Moreover, hyperoxic BOLD modulation is primarily sensitive to CBV_v,0 and largely insensitive to OEF₀; in fact, hyperoxia can be used to estimate CBV_v,0.⁶² The oxygen saturation change due to hyperoxia stimulus is independent of the baseline oxygen saturation over most of the physiological range. In contrast, the oxygen saturation change to a hypercapnic challenge is linearly related to the resting saturation. The sensitivity pattern of the hyperoxic modulation makes the estimation of OEF₀ with the new framework completely reliant on the flow-diffusion model approximations that link CBV_v,0 to OEF₀. The model approximations are indeed less influential with hypercapnia, which has a larger sensitivity to OEF₀ with respect to CBV_v,0, making the hypercapnia approach less noisy and biased.

Limitations of the method

The main limitations of the new method are mostly shared with dc-fMRI.¹⁴ The approach using hypercapnia relies on a local CBF increase, a vascular reserve, which may be absent in diseases such ischemic stroke, where vessels may be maximally dilated in an attempt to maintain perfusion. In addition, the method might be vulnerable to larger than expected changes in ρ or k, which are probably not independent. Although large changes in these parameters appear unlikely in many brain diseases, we might expect relevant tissue and vascular remodeling in some disease, such as brain tumors.⁶³ The limitations of the approach in diseases with concurrent very high P_mO_2,0 and long MCTT₀ are noted earlier and should be assessed in future studies.

Study limitations

The main limitation of the simulation study lies in the assumption of an exact analytical model with the error in the estimate of OEF₀ being introduced only by the limited number of measurable variables. The main simplifying assumption of the model was the replacement of the CTT in one straight capillary with the MCTT in the voxel capillary bed. This is an approximation, since, due to the non-linear mapping between CTT, OEF and oxygen diffusion between capillary and tissue, the complete CTT distribution within a capillary bed affects the macroscopic OEF.⁶⁴ Without changing MCTT, OEF can be increased through homogenization of the CTT among capillaries. Future extension of the model might include the CTT heterogeneity (CTTH), a measure of the second moment of the CTT distribution within the capillary bed.⁶⁵

With respect to the in-vivo validation using dc-fMRI, a limiting factor was related to the investigation of the proposed model proportionality constant A · ρ/k and P_mO_2,0. These estimates were evaluated assuming the OEF₀ derived from the dc-fMRI machine learning analysis²² to be exact. In fact, noise in the dc-fMRI OEF₀ limited our investigation of the model parameters to global evaluation within the GM. Moreover, another limitation was the problem of having two unknowns and one equation. By exploiting the different effects of these parameters on the non-linear mapping between variables, we were able to get insight into both parameters, however only at a between subjects’ average level. Alternative approaches should be used to investigate the different physiological parameters (e.g., ρ and k) contributing to the proportionality constant, which cannot be separately investigated using standard fMRI approaches. Comparison against non-MRI technology would be essential for definitive validation of the approach. Future validation is also necessary beyond the healthy controls involved in the study, to populations affected by diseases that might alter brain metabolism and for which the proposed model’s validity might reduce.

Conclusion

We introduced a novel single gas calibrated fMRI framework integrating a steady-state flow-diffusion model of oxygen transport into the BOLD signal model. Uncertainty in the integrated model is driven by variability in a proportionality constant A · ρ/k, that depends on tissue and microvascular structure at a fixed field strength, and PmO_2,0. The advantage of the new framework lies in the limited influence of these parameters on OEF₀. Even by fixing them to plausible values, the simulations showed the OEF₀ RMSE was below 0.05 over a wide range of physiology meaning that the method may reliably identify between-subjects OEF₀ differences greater than approximately 10%. Only with concurrently very high PmO_2,0 (>25 mmHg) and long MCTT₀ (>2.5 s) did the RMSE increase significantly suggesting that the method should work in diseases not drastically altering brain physiology. Importantly, the approach was highly robust to measurement noise. The method, when using hypercapnia, compared well with dc-fMRI and with whole-brain OEF₀ derived using the OxFlow method. Lack of positive results when using hyperoxia may be related to the method being primarily sensitive to CBV_v,0. The simplified calibrated fMRI method using hypercapnia has potential for application in clinical settings.

Table 1.

Main variables and abbreviations used in the study, reported in alphabetical order.

Abbreviation	Meaning	Units	Abbreviation	Meaning	Units
₀	As subscript defines the physiological variable at baseline	/	[Hb]	Concentration of hemoglobin in blood	g/dL
α	Grubb exponent relating fractional change in CBV_v to fractional change in CBF	Dimensionless	β	Field strength and vessel geometry dependent exponent within the steady-state BOLD signal model	Dimensionless
β-weight	Coefficient of the GLM	/	γ	Gyromagnetic ratio of the proton	rad/s/T
Δχ_do	Magnetic susceptibility difference between fully oxygenated and fully deoxygenated blood	Relative	ΔBOLD/BOLD₀	Relative change in BOLD signal	Relative
ΔS	Tag-Control ASL image	Not defined	ε	Oxygen plasma solubility	mL/mmHg/dL
η	Tagging inversion efficiency of PCASL	Dimensionless	η _INV	Scaling factor accounting for reduction in tagging efficiency due to background suppression	Dimensionless
λ	Water partition coefficient of the tissue	mL/g	ρ	CBV_v,0/CBV_cap,0	Relative
φ	Oxygen binding capacity of hemoglobin	mL/g	Φ	Phase MRI image	rad
τ	Labelling duration of PCASL	s	A	Field strength and vessel geometry dependent proportionality constant within the steady-state BOLD signal model	s⁻¹g^−βdL^β
ASL	Arterial spin labelling	/	B₀	Static magnetic field	T
BOLD	Blood oxygen level dependent	Not Defined	C_aO₂	Concentration of oxygen in arteries	mL/dL
CBF	Cerebral blood flow	mL/100 g/min	CBF/CBF₀	Fractional change in CBF	Relative
CBV_cap	Capillary blood volume	Relative (or ml/100 g)	CBV_v	dHb-sensitive blood volume	Relative (or ml/100 g)
CMRO₂	Cerebral metabolic rate of oxygen	μmol/100 g/min	CBV_v/CBV_v,0	Fractional change in CBV_v	Relative
CSF	Cerebrospinal fluid	/	CTT	Capillary transit time	s
CoV	Coefficient of variation	Relative	CVR	Cerebrovascular reactivity ( CBF/CBF₀/mmHg of CO₂)	%/mmHg
dc-fMRI	Dual-Calibrated functional MRI	/	dHb	Deoxy-hemoglobin concentration in tissue	g/100 g
GLM	General Linear Model	/	GM	Grey matter	/
GRE	Gradient Echo Sequence	/	k	Effective permeability to oxygen of the capillary endothelium and brain tissue	μmol/mmHg/mL/min
h	Hill constant involved in the non-linear relationship between oxygen partial pressure and hemoglobin saturation in blood	Dimensionless	hc-fMRI+	Single gas calibrated fMRI using a hypercapnic modulation and the steady-state BOLD signal model extended with the proposed flow-diffusion analytical framework of oxygen transport	/
hc-fMRI	Single gas calibrated fMRI using a hypercapnic modulation and the steady-state BOLD signal model	/	Hct	Hematocrit	%
ho-fMRI+	Single gas calibrated fMRI using a hyperoxic modulation and the steady-state BOLD signal model extended with the proposed flow-diffusion analytical framework of oxygen transport		ho-fMRI	Single gas calibrated fMRI using a hyperoxic modulation and the steady-state BOLD signal model	/
MCTT	Mean capillary transit time	s	M	Maximum BOLD modulation	Relative
OxFlow	Validated macrovascular global measure of OEF₀, inferred through phase measures of the magnetic susceptibility of blood in the sagittal sinus relative to surrounding tissue	/	OEF	Oxygen Extraction Fraction	Relative
P₅₀	Oxygen partial pressure when half of hemoglobin is saturated with oxygen	mmHg	P_aO₂	Partial pressure of oxygen in arteries	mmHg
P_aCO₂	Partial pressure of carbon dioxide in arteries	mmHg	P_capO₂	Partial pressure of oxygen in the capillary	mmHg
PCASL	Pseudo-continuous ASL	/	P_ETO₂	End-tidal partial pressure of oxygen	mmHg
P_ETCO₂	End-tidal partial pressure of carbon dioxide	mmHg	P_mO₂	Partial pressure of oxygen at the mitochondria	mmHg
PLD	Post label delay of PCASL	s	R₂*\|_dHb	Rate of free induction decay due to dHb	1/s
RMSE	Root mean square error	Variable	S₀	Proton Density Image for ASL normalization	Not defined
S_capO₂	Capillary oxygen saturation of hemoglobin	Relative	S_aO₂	Arterial oxygen saturation of hemoglobin	Relative
SD	Standard deviation	Variable	SE	Standard error (standard deviation of the mean)	Variable
S_vO₂	Venous oxygen saturation of hemoglobin	Relative	SNR	Signal to Noise Ratio	Dimensionless
SSS	Superior Sagittal Sinus	/	T1_b	MRI longitudinal relaxation time constant of blood	s
TE	Time of echo of the MRI sequence	s	TOF	Time of flight MRI	/
TR	Time of repetition of the MRI sequence		WM	White matter	/

Supplemental Material

sj-pdf-1-jcb-10.1177_0271678X221077332 - Supplemental material for A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO₂) with single gas calibrated fMRI

Supplemental material, sj-pdf-1-jcb-10.1177_0271678X221077332 for A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO₂) with single gas calibrated fMRI by Antonio M Chiarelli, Michael Germuska, Hannah Chandler, Rachael Stickland, Eleonora Patitucci, Emma Biondetti, Daniele Mascali, Neeraj Saxena, Sharmila Khot, Jessica Steventon, Catherine Foster, Ana E Rodríguez-Soto, Erin Englund, Kevin Murphy, Valentina Tomassini, Felix W Wehrli and Richard G Wise in Journal of Cerebral Blood Flow & Metabolism

Footnotes

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: MG thanks the Wellcome Trust for its support via a Sir Henry Dale Fellowship (220575/Z/20/Z). This project was partially supported by the UK Engineering and Physical Sciences Research Council (EP/S025901/1). This work was supported by the National Institute of Health (P41-EB029460). HLC and MG were funded by a Wellcome Strategic Award to CUBRIC, ‘Multi-scale and multi-modal assessment of coupling in the healthy and diseased brain’ (104943/Z/14/Z). RCS, EP and CF were supported by Wellcome PhD studentships. KM was supported by a Wellcome Senior Research Fellowship, “Assessing the health of ageing blood vessels in the brain using fMRI” (200804/Z/16/Z). The study was partially supported by the MS Society UK.

Acknowledgments

This work was partially conducted under the framework of the Departments of Excellence 2018–2022 initiative of the Italian Ministry of Education, University and Research for the Department of Neuroscience, Imaging and Clinical Sciences (DNISC) of the University of Chieti-Pescara, Italy. For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission. The data that support the findings of this study are available upon reasonable request.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Authors’ contributions

AMC and MG developed imaging and analysis methods, and analyzed and interpreted the data. AMC, MG and RGW drafted the manuscript. FWW, RGW, KM and VT conceived the project. MG, HLC, RCS, EP, NS, SK, SJ, CF and KM set up and executed the experiment. DM, EB, AERS and EE contributed to data analysis. All authors revised and approved the final submission.

ORCID iDs

Antonio M Chiarelli

Michael Germuska

Emma Biondetti

Supplemental material

Supplemental material for this article is available online.

References

Magistretti PJ and Allaman I. Brain Energy Metabolism. In: Pfaff D.W. (eds) Neuroscience in the 21st Century From Basic to Clinical. 1st ed. New York: Springer, 2013, pp. 1591–1620.

Barhoum

Langham

Magland

, et al. Method for rapid MRI quantification of global cerebral metabolic rate of oxygen. J Cereb Blood Flow Metab 2015; 35: 1616–1622.

Bolar

Rosen

Sorensen

, et al. QUantitative imaging of eXtraction of oxygen and TIssue consumption (QUIXOTIC) using venular-targeted velocity-selective spin labeling. Magn Reson Med 2011; 66: 1550–1562.

Bulte

Kelly

Germuska

, et al. Quantitative measurement of cerebral physiology using respiratory-calibrated MRI. NeuroImage 2012; 60: 582–591.

Fan

Govindarajan

Kinkel

, et al. Quantitative oxygen extraction fraction from 7-Tesla MRI phase: Reproducibility and application in multiple sclerosis. J Cereb Blood Flow Metab 2015; 35: 131–139.

Guo

Wong

EC.

Venous oxygenation mapping using velocity-selective excitation and arterial nulling. Magn Reson Med 2012; 68: 1458–1471.

Yablonskiy

DA.

Quantitative BOLD: Mapping of human cerebral deoxygenated blood volume and oxygen extraction fraction: default state. Magn Reson Med 2007; 57: 115–126.

Grgac

, et al. Calibration and validation of TRUST MRI for the estimation of cerebral blood oxygenation. Magn Reson Med 2012; 67: 42–49.

Barhoum

Rodgers

Langham

, et al. Comparison of MRI methods for measuring whole-brain venous oxygen saturation. Magn Reson Med 2015; 73: 2122–2128.

10.

Zhu

Yablonskiy

DA.

Validation of oxygen extraction fraction measurement by qBOLD technique. Magn Reson Med 2008; 60: 882–888.

11.

Zhang

Liu

Gupta

, et al. Quantitative mapping of cerebral metabolic rate of oxygen (CMRO2) using quantitative susceptibility mapping (QSM). Magn Reson Med 2015; 74: 945–952.

12.

Gauthier

Hoge

RD.

Magnetic resonance imaging of resting OEF and CMRO₂ using a generalized calibration model for hypercapnia and hyperoxia. NeuroImage 2012; 60: 1212–1225.

13.

Germuska

Chandler

Stickland

, et al. Dual-calibrated fMRI measurement of absolute cerebral metabolic rate of oxygen consumption and effective oxygen diffusivity. NeuroImage 2019; 184: 717–728.

14.

Germuska

Wise

RG.

Calibrated fMRI for mapping absolute CMRO2: Practicalities and prospects. NeuroImage 2019; 187: 145–153.

15.

Wise

Harris

Stone

, et al. Measurement of OEF and absolute CMRO2: MRI-based methods using interleaved and combined hypercapnia and hyperoxia. NeuroImage 2013; 83: 135–147.

16.

Davis

Kwong

Weisskoff

, et al. Calibrated functional MRI: mapping the dynamics of oxidative metabolism. Proc Natl Acad Sci U S A 1998; 95: 1834–1839.

17.

Champagne

Coverdale

Germuska

, et al. Changes in volumetric and metabolic parameters relate to differences in exposure to Sub-concussive head impacts. J Cereb Blood Flow Metab 2020; 40: 1453–1467.

18.

De Vis

Petersen

Bhogal

, et al. Calibrated MRI to evaluate cerebral hemodynamics in patients with an internal carotid artery occlusion. J Cereb Blood Flow Metab 2015; 35: 1015–1023.

19.

Lajoie

Nugent

Debacker

, et al. Application of calibrated fMRI in Alzheimer’s disease. NeuroImage Clin 2017; 15: 348–358.

20.

Merola

Germuska

Murphy

, et al. Assessing the repeatability of absolute CMRO2, OEF and haemodynamic measurements from calibrated fMRI. NeuroImage 2018; 173: 113–126.

21.

Merola

Germuska

Warnert

, et al. Mapping the pharmacological modulation of brain oxygen metabolism: the effects of caffeine on absolute CMRO2 measured using dual calibrated fMRI. NeuroImage 2017; 155: 331–343.

22.

Germuska

Chandler

Okell

, et al. A Frequency-Domain machine learning method for dual-calibrated fMRI mapping of oxygen extraction fraction (OEF) and cerebral metabolic rate of oxygen consumption (CMRO2). Front Artif Intell 2020; 3. Epub ahead of print DOI: 10.3389/frai.2020.00012.

23.

Gjedde

Cerebral blood flow change in arterial hypoxemia is consistent with negligible oxygen tension in brain mitochondria. NeuroImage 2002; 17: 1876–1881.

24.

Hayashi

Watabe

Kudomi

, et al. A theoretical model of oxygen delivery and metabolism for physiologic interpretation of quantitative cerebral blood flow and metabolic rate of oxygen. J Cereb Blood Flow Metab off J Tab 2003; 23: 1314–1323.

25.

Lee

Langham

Rodriguez-Soto

, et al. Multiplexed MRI methods for rapid estimation of global cerebral metabolic rate of oxygen consumption. NeuroImage 2017; 149: 393–403.

26.

Boxerman

Bandettini

Kwong

, et al. The intravascular contribution to FMRI signal change: Monte Carlo modeling and diffusion-weighted studies in vivo. Magn Reson Med 1995; 34: 4–10.

27.

Hoge

Atkinson

Gill

, et al. Investigation of BOLD signal dependence on cerebral blood flow and oxygen consumption: the deoxyhemoglobin dilution model. Magn Reson Med 1999; 42: 849–863.

28.

Buxton

RB.

Introduction to functional magnetic resonance imaging: Principles and techniques. Cambridge: Cambridge University Press, 2009.

29.

Germuska

Merola

Murphy

, et al. A forward modelling approach for the estimation of oxygen extraction fraction by calibrated fMRI. NeuroImage 2016; 139: 313–323.

30.

Hyder

Shulman

Rothman

DL.

A model for the regulation of cerebral oxygen delivery. J Appl Physiol (1985) 1998; 85: 554–564.

31.

Ogawa

Menon

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 1993; 64: 803–812.

32.

Blockley

Griffeth

VEM

Buxton

RB.

A general analysis of calibrated BOLD methodology for measuring CMRO2 responses: comparison of a new approach with existing methods. NeuroImage 2012; 60: 279–289.

33.

Blockley

Griffeth

VEM

Simon

, et al. Calibrating the BOLD response without administering gases: Comparison of hypercapnia calibration with calibration using an asymmetric spin echo. NeuroImage 2015; 104: 423–429.

34.

Blockley

Stone

AJ.

Improving the specificity of R2′ to the deoxyhaemoglobin content of brain tissue: Prospective correction of macroscopic magnetic field gradients. NeuroImage 2016; 135: 253–260.

35.

Hua

Liu

Kim

, et al. MRI techniques to measure arterial and venous cerebral blood volume. NeuroImage 2019; 187: 17–31.

36.

Gjedde

Poulsen

Østergaard

On the oxygenation of hemoglobin in the human brain. In: Eke A, Delpy DT (eds) Oxygen transport to tissue XXI. Boston, MA: Springer US, pp. 67–81.

37.

Gjedde

Johannsen

Cold

, et al. Cerebral metabolic response to low blood flow: possible role of cytochrome oxidase inhibition. J Cereb Blood Flow Metab 2005; 25: 1183–1196.

38.

Okell

Chappell

Kelly

, et al. Cerebral blood flow quantification using vessel-encoded arterial spin labeling. J Cereb Blood Flow Metab 2013; 33: 1716–1724.

39.

Schmithorst

Hernandez-Garcia

Vannest

, et al. Optimized simultaneous ASL and BOLD functional imaging of the whole brain. J Magn Reson Imaging JMRI 2014; 39: 1104–1117.

40.

Jain

Langham

Wehrli

FW.

MRI estimation of global brain oxygen consumption rate. J Cereb Blood Flow Metab 2010; 30: 1598–1607.

41.

Langham

Magland

Epstein

, et al. Accuracy and precision of MR blood oximetry based on the long paramagnetic cylinder approximation of large vessels. Magn Reson Med 2009; 62: 333–340.

42.

Langham

Magland

Floyd

, et al. Retrospective correction for induced magnetic field inhomogeneity in measurements of large-vessel hemoglobin oxygen saturation by MR susceptometry. Magn Reson Med 2009; 61: 626–633.

43.

Gai

Taki

Kato

, et al. Regulation of hemoglobin affinity for oxygen by carbonic anhydrase. J Lab Clin Med 2003; 142: 414–420.

44.

Jenkinson

Beckmann

Behrens

TEJ

, et al. FSL. NeuroImage 2012; 62: 782–790.

45.

Jenkinson

Bannister

Brady

, et al. Improved optimization for the robust and accurate linear registration and motion correction of brain images. NeuroImage 2002; 17: 825–841.

46.

Smith

SM.

Fast robust automated brain extraction. Hum Brain Mapp 2002; 17: 143–155.

47.

Zhang

Brady

Smith

Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans Med Imaging 2001; 20: 45–57.

48.

Alsop

Detre

Golay

, et al. Recommended implementation of arterial spin-labeled perfusion MRI for clinical applications: a consensus of the ISMRM perfusion study group and the European consortium for ASL in dementia. Magn Reson Med 2015; 73: 102–116.

49.

Mutsaerts

HJMM

Steketee

RME

Heijtel

DFR

, et al. Inter-vendor reproducibility of pseudo-continuous arterial spin labeling at 3 tesla. Plos One 2014; 9: e104108.

50.

Pilkinton

Hiraki

Detre

, et al. Absolute cerebral blood flow quantification with pulsed arterial spin labeling during hyperoxia corrected with the simultaneous measurement of the longitudinal relaxation time of arterial blood. Magn Reson Med 2012; 67: 1556–1565.

51.

Friston

KJ.

Statistical parametric mapping. In: Functional neuroimaging: Technical foundations. San Diego, CA: Academic Press, 1994, pp. 79–93.

52.

Insiripong

Supattarobol

Jetsrisuparb

Comparison of hematocrit/hemoglobin ratios in subjects with alpha-thalassemia, with subjects having chronic kidney disease and normal subjects. Southeast Asian J Trop Med Public Health 2013; 44: 707–711.

53.

Spees

Yablonskiy

Oswood

, et al. Water proton MR properties of human blood at 1.5 tesla: magnetic susceptibility, T1, T2, T, and non-Lorentzian signal behavior. Magn Reson Med 2001; 45: 533–542.

54.

Kastrup

Krüger

Neumann-Haefelin

, et al. Assessment of cerebrovascular reactivity with functional magnetic resonance imaging: comparison of CO₂ and breath holding. Magn Reson Imaging 2001; 19: 13–20.

55.

Germuska

Stickland

Chiarelli

, et al. Quantitative mapping of cerebral oxygen metabolism using breath-hold calibrated fMRI. bioRxiv 2021; 2021.04.08.438939.

56.

Wilson

DF.

Regulation of cellular metabolism: programming and maintaining metabolic homeostasis. J Appl Physiol (1985) 2013; 115: 1583–1588.

57.

Buxton

RB.

The thermodynamics of thinking: connections between neural activity, energy metabolism and blood flow. Philos Trans R Soc Lond B Biol Sci 2021; 376: 20190624.

58.

Driver

Wharton

Croal

, et al. Global intravascular and local hyperoxia contrast phase-based blood oxygenation measurements. NeuroImage 2014; 101: 458–465.

59.

Liu

Test-retest reproducibility of a rapid method to measure brain oxygen metabolism. Magn Reson Med 2013; 69: 675–681.

60.

Jiang

Deng

Franklin

, et al. Validation of T2-based oxygen extraction fraction measurement with 15O positron emission tomography. Magn Reson Med 2021; 85: 290–297.

61.

Perlmutter

Powers

Herscovitch

, et al. Regional asymmetries of cerebral blood flow, blood volume, and oxygen utilization and extraction in normal subjects. J Cereb Blood Flow Metab 1987; 7: 64–67.

62.

Blockley

Griffeth

VEM

Germuska

, et al. An analysis of the use of hyperoxia for measuring venous cerebral blood volume: Comparison of the existing method with a new analysis approach. NeuroImage 2013; 72: 33–40.

63.

Owen

Alarcón

Maini

, et al. Angiogenesis and vascular remodelling in normal and cancerous tissues. J Math Biol 2009; 58: 689–721.

64.

Wei

Wang

RK.

Capillary flow homogenization during functional activation revealed by optical coherence tomography angiography based capillary velocimetry. Sci Rep 2018; 8: 4107.

65.

Jespersen

Østergaard

The roles of cerebral blood flow, capillary transit time heterogeneity, and oxygen tension in brain oxygenation and metabolism. J Cereb Blood Flow Metab 2012; 32: 264–277.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.36 MB

A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO 2 ) with single gas calibrated fMRI

Abstract

Keywords

Introduction

Methods

Analytical modeling

BOLD model and the Dual-Calibrated fMRI experiment

Flow-Diffusion model of oxygen transport

Integration of the Flow-Diffusion model of oxygen transport into the BOLD model for calibrated-fMRI quantification of CMRO2

Simulations

MRI experiment

fMRI data processing

Gas recordings processing

fMRI processing

Dual-calibrated fMRI analysis

Single gas calibrated fMRI analysis

OxFlow data processing

Statistical analysis

Results

Simulations

In-Vivo evaluation of gas, CBF and BOLD modulations

In-Vivo estimation of modeling parameters

In-vivo estimation of oxygen extraction fraction: Calibrated fMRI vs. dual-calibrated fMRI (section 3)

In-vivo estimation of oxygen extraction fraction: Calibrated fMRI vs. OxFlow (section 4)

Discussion

Simulations

Modeling parameters

Comparison with dc-fMRI and OxFlow

Limitations of the method

Study limitations

Conclusion

Supplemental Material

sj-pdf-1-jcb-10.1177_0271678X221077332 - Supplemental material for A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO2) with single gas calibrated fMRI

Footnotes

Funding

Acknowledgments

Declaration of conflicting interests

Authors’ contributions

ORCID iDs

Supplemental material

References

Supplementary Material

Integration of the Flow-Diffusion model of oxygen transport into the BOLD model for calibrated-fMRI quantification of CMRO₂

sj-pdf-1-jcb-10.1177_0271678X221077332 - Supplemental material for A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO₂) with single gas calibrated fMRI