Quantitative functional BOLD (qfBOLD): A combined gradient-echo and spin-echo framework for oxygen extraction fraction (OEF) mapping with functional MRI

Abstract

We present a novel framework for OEF mapping with MRI, based on temporal variations in gradient echo (GE) and spin-echo (SE) BOLD signals induced by isometabolic modulations in CBF. This approach, termed quantitative functional BOLD (qfBOLD), exploits dynamic variations in relaxation times rather than measuring baseline values as in qBOLD, thereby isolating deoxyhaemoglobin (dHb) effects. The interaction between dHb-induced extravascular field distortions and water diffusion allows for decoupling OEF and dHb-sensitive cerebral blood volume with a single modulation in brain physiology. Furthermore, the method avoids functional CBF measures via arterial spin labelling which is required by calibrated (c)fMRI. This advancement may enhance signal-to-noise ratio and spatiotemporal resolution, making qfBOLD applicable to both grey matter (GM) and white matter (WM). Monte Carlo simulations were used to investigate the method. In vivo feasibility assessment using a hypercapnic breath-holding task yielded OEF values of 37.0% ± 2.9% and 41.6% ± 2.9% in GM and WM, respectively, and significant correlations with cfMRI in GM (qfBOLD vs cfMRI r = 0.71, p < 10⁻³) and with relaxometry-based measures in the superior sagittal sinus (GM qfBOLD vs TRUST r = 0.51, p < 0.05, WM qfBOLD vs TRUST r = 0.61, p < 0.01). Future efforts will aim to improve the method’s accuracy by attenuating intravascular signals and by refining WM modelling.

Keywords

Calibrated functional (f)MRI gradient-echo blood oxygen level dependent (BOLD) fMRI hypercapnia oxygen extraction fraction (OEF)spin-echo BOLD fmri

Introduction

Oxidative metabolism provides most of the brain’s energy.^1,2 The brain lacks energetic reserves, and it depends on oxygen supply via cerebral blood flow (CBF).^2,3 The brain’s use of oxygen is quantified by the cerebral metabolic rate of oxygen (CMRO₂).^4–6 Oxygen transfer must obey the principle of conservation of mass which states that CMRO₂ equals the rate of exit of oxygen from the vascular compartment.² This implies that CMRO₂ is equal to the product of oxygen supply (itself product of CBF and arterial oxygen concentration (CaO₂)) and the fraction of oxygen extracted (i.e. the oxygen extraction fraction (OEF)).⁷

Magnetic resonance imaging (MRI) can be used to quantify CMRO₂, offering advantages over the gold-standard approach of ¹⁵O PET, including the avoidance of ionising radiation and the need for a cyclotron to make short-lived contrast agents.^8,9 While different MRI techniques can quantify perfusion (e.g. arterial spin labelling (ASL)),¹⁰ measuring OEF is more challenging. It requires the estimation of deoxyhaemoglobin (dHb) concentration in blood ([dHb]) to infer venous saturation (SvO₂) and OEF.¹¹ OEF can be reliably measured in large veins, but such methods do not localise the site of oxygen extraction.^12,13

In contrast, evaluation of [dHb] in microvasculature allows the generation of OEF maps but it is affected by two main limitations. Firstly, estimating [dHb] requires knowledge of the dHb-sensitive cerebral blood volume (CBV_dHb).¹⁴ Secondly, measurements can be affected by non-blood susceptibility sources, such as iron and myelin. This is mainly a concern for methods such as quantitative blood oxygen level dependent (qBOLD) MRI and quantitative susceptibility mapping (QSM), which aim to quantify OEF from an evaluation of signal decay and phase changes with echo time (TE), respectively.^15–18 Calibrated functional (cf)MRI addresses the latter problem by probing the temporal modulations in the transverse relaxation rate (R₂*, relaxation time T₂* = 1/R₂*) related to dHb variations. It combines gradient echo BOLD signal changes (BOLD_GE) and ASL functional recordings with isometabolic modulations of CBF to infer, at a given TE, the maximum BOLD_GE increase obtainable with complete removal of dHb.¹⁹ The most informative approach is dual cfMRI, which uses a complex protocol alternating between hypercapnia (a rise in CO₂ in arterial blood increasing CBF) and hyperoxia (increased O₂ in arterial blood) through gas inhalation, to decouple baseline OEF from CBV_dHb and thereby map GM CMRO₂.^20–24 The principal shortcomings of cfMRI methods are the practical complexity and the need to correlate BOLD and ASL perfusion modulations. ASL has low temporal resolution and temporal signal-to-noise ratio (tSNR). In practice, it is generally feasible only in GM and not in WM. These features limit the applicability of cfMRI methods.²¹

Here we present a novel fMRI framework that estimates baseline [dHb] and OEF by combining T₂*-weighted BOLD_GE and T₂-weighted spin echo BOLD (BOLD_SE) fMRI concurrently acquired during modulation of CBF. The method estimates OEF with the selective sensitivity to dHb of cfMRI but with a simplified protocol that requires hypercapnia only, and without the need for concurrent ASL recordings. The absence of functional ASL should produce measurements with higher SNR with the possibility to increase spatiotemporal resolution. Moreover, the approach may also be applied to WM where ASL is challenging. We refer to the method as quantitative functional BOLD (qfBOLD), since it uses functional BOLD_SE and BOLD_GE signal modulations to infer relaxation-rate changes and quantify OEF, analogous to qBOLD in its use of SE and GE deoxyhaemoglobin-sensitive signal modelling. However, unlike qBOLD, which typically relies on baseline multi-TE signal behaviour to estimate absolute relaxation properties such as T₂*, T₂ and T₂′ (1/T₂′ = 1/T₂* − 1/T₂), qfBOLD exploits temporal fluctuations in BOLD_SE and BOLD_GE signals assessed at a single TE.

Here, we firstly provide a descriptive overview of the approach. Secondly, we introduce the qfBOLD analytical framework, where the simple but robust and interpretable Davis model of the BOLD signal is used.¹⁹ Thirdly, we report on the validation of the method based on Monte Carlo simulations involving a multi-compartmental model of the BOLD_SE and BOLD_GE signals, exploring effects related to vessel topology, contribution of the intravascular signal and experimental noise. Finally, we test the method in healthy human subjects. We implemented a tailored BOLD–ASL sequence to concurrently acquire BOLD_GE, BOLD_SE and pseudo-continuous ASL (pCASL) data during a breath-holding task (breath-hold, BH) that induced CBF changes via hypercapnic modulation. The sequence was developed specifically to compare the qfBOLD method to a single cfMRI approach that we recently developed to map OEF.^7,25 qfBOLD was also compared to independent macrovascular measures of OEF derived from the superior sagittal sinus (SSS, TRUST, T₂ relaxation under spin tagging).¹²

Methods

Overview of qfBOLD

The qfBOLD method estimates OEF by combining T₂*-weighted BOLD_GE and T₂-weighted spin echo BOLD_SE fMRI acquired during modulation of CBF. The framework exploits the fact that, in the presence of microvasculature and water diffusion, both extravascular R₂ (1/T₂) and R₂* (1/T₂*) are linearly related to CBV_dHb but they scale non-linearly, and differently, with respect to [dHb].^26–28 The effect is stronger for smaller vascular compartments, such as capillaries. Since the BOLD signal modulations are approximately proportional to relaxation rates changes (−ΔR₂* for BOLD_SE and −ΔR₂ for BOLD_SE), they will be non-linear functions of changes in [dHb] (Δ[dHb]) and will instead preserve a linear behaviour with CBV_dHb and its changes. Crucially, since the BOLD_SE weighting has an increased sensitivity to diffusion effects compared to BOLD_GE and it primarily ‘senses’ the capillary compartment, the BOLD_SE signal change exhibits an accentuated supralinear dependence on ∆[dHb] compared to BOLD_GE. This implies that the contribution of ∆[dHb] to the signals can be isolated by comparing the BOLD_GE and the BOLD_SE signal changes, specifically, by taking their ratio. Notably, the functions linking BOLD_SE and BOLD_GE modulations to Δ[dHb] are influenced by vessel topology (for a fixed echo time, TE, readout scheme and intravascular signal suppression).^29–31 The sensitivity of the BOLD_SE/BOLD_GE signal change ratio to vessel topology is particularly strong when large vessels (>100 µm) are considered.³² However, we speculate that the effect of ∆[dHb] on the BOLD_SE/BOLD_GE signal change ratio is generally stronger than the effect induced by variability in the average vessel sizes within a voxel where only microvasculature is present. This is because capillary and venular topology as well as their volume ratio within each voxel are not expected to vary extensively.^33–35

Importantly, analogously to cfMRI relying on hypercapnia, a modulation in CBF makes ∆[dHb] a function of baseline [dHb],²³ implying that the comparison of the BOLD_SE and BOLD_GE modulations during vasodilation is sensitive to baseline venous saturation and OEF.²³ When other mechanisms are used to alter vascular susceptibility, like gadolinium administration or hyperoxia, the method does not depend on baseline [dHb], thus becoming selectively sensitive to vessel topology (i.e. vessel size imaging).^36–39 When using CBF modulations to introduce sensitivity to baseline [dHb], the CBF changes affect the BOLD signals in a non-linear manner. However, the non-linear contributions of CBF modulations differ only slightly between BOLD_GE and BOLD_SE signal changes. As a result, the ratio of BOLD_SE to BOLD_GE modulations is not only independent of CBV_dHb, but also nearly independent of the extent of CBF modulation, removing the need for ASL measurement of CBF.

qfBOLD analytical modelling

The rate of signal decay of BOLD_GE due to dHb within a voxel with microvasculature can be expressed, according to the Davis model, as¹⁹:

R_{2}^{*} |_{d H b} = A_{G E} \cdot C B V_{d H b} \cdot {[d H b]}^{β_{G E}} .

(1)

Here, $β$ _GE encodes the supralinear dependence of R₂* on [dHb] and it depends (as A_GE) on the MRI sequence (TE, readout scheme), water diffusion around microvessels, vascular topology (e.g. average vessel radius) and intravascular/extravascular contributions to the BOLD_GE signal. Values of β_GE ≈ 1.3–1.5 have been reported in the literature for BOLD_GE at 3 T with TE_GE ≈ 30 ms.^7,19 The net supralinear behaviour primarily arises from a linear large vessel component combined with supralinear contributions from capillaries to the extravascular relaxivity.^7,19

In the presence of small ∆ $R_{2}^{*}$ induced by dHb modulation, the BOLD_GE signal change can be expressed as^19,22:

\begin{array}{l} \frac{Δ B O L D_{G E}}{B O L D_{G E}} = - T E_{G E} \cdot Δ R_{2}^{*} = \\ M_{G E} \cdot {1 - (\frac{C B V_{dHb, m}}{C B V_{d H b}}) \cdot {(\frac{[{dHb}_{m}]}{[dHb]})}^{β_{G E}}} \end{array}

(2)

where the subscript m indicates a modulation relative to the baseline value and M_GE is the maximum BOLD_GE signal change that can be obtained with complete removal of dHb. M_GE depends on baseline quantities and is equal to:

M_{G E} = T E_{G E} \cdot A_{G E} \cdot C B V_{d H b} \cdot {([d H b])}^{β_{G E}}

(3)

$C B V_{dHb}$ changes can be expressed as a function of CBF modulations via the Grubb relation,^41,42 whereas [dHb] changes can be expressed as a function of CBF and CMRO₂ changes obtaining:

\frac{Δ B O L D_{G E}}{B O L D_{G E}} = M_{G E} \cdot {1 - {(\frac{{CBF}_{m}}{CBF})}^{α - β_{G E}} \cdot {(\frac{{CMRO}_{2, m}}{{CMRO}_{2}})}^{β_{G E}}}

(4)

with α the Grubb exponent (α = 0.18).⁴¹

Assuming isometabolic vasodilation, equation (4) becomes:

\begin{array}{l} \frac{Δ B O L D_{G E}}{B O L D_{G E}} = M_{G E} \cdot {1 - {(\frac{{CBF}_{m}}{CBF})}^{α - β_{G E}}} = \\ M_{G E} \cdot {1 - {(1 + \frac{Δ CBF}{CBF})}^{α - β_{G E}}} \end{array}

(5)

If flow changes are of limited size, equation (5) can be linearised thereby obtaining⁴²:

\frac{Δ B O L D_{G E}}{B O L D_{G E}} = M_{G E} \cdot (β_{G E} - α) \cdot \frac{Δ C B F}{C B F}

(6)

The refocussing pulse in a BOLD_SE sequence reduces the effect of dHb on the BOLD signal while enhancing BOLD signal sensitivity to water diffusion and capillary extravascular effects.⁴³ Within the Davis model framework, the refocussing effect may be conceived as lowering the multiplicative constant A and increasing the exponent factor β. Consequently, we can derive a comparable set of equations for BOLD_SE related to small perturbations of R₂ due to flow changes, using different values of A and β compared to BOLD_GE, resulting in:

\frac{Δ B O L D_{S E}}{B O L D_{S E}} = M_{S E} \cdot (β_{S E} - α) \cdot \frac{Δ C B F}{C B F}

(7)

With M_SE equal to:

M_{S E} = T E_{S E} \cdot A_{S E} \cdot C B V_{d H b} \cdot {([d H b])}^{β_{S E}}

(8)

By comparing BOLD_SE to BOLD_GE signal changes we obtain:

\begin{array}{l} \frac{\frac{Δ B O L D_{S E}}{B O L D_{S E}}}{\frac{Δ B O L D_{G E}}{B O L D_{G E}}} = \frac{T E_{S E} \cdot A_{S E} \cdot (β_{S E} - α)}{T E_{G E} \cdot A_{G E} \cdot (β_{G E} - α)} \cdot \\ {([d H b])}^{β_{S E - G E}} = C \cdot {([d H b])}^{β_{S E - G E}} \end{array}

(9)

where C is a multiplicative constant and $β_{S E - G E}$ (value larger than 0) is the exponent representing the difference in $β s$ between BOLD_SE and BOLD_GE.

Hence, the BOLD_SE/BOLD_GE signal change ratio can be assumed to be independent of baseline CBV_dHb and CBF modulations and, importantly, is a monotonic function of [dHb]. SvO₂ and OEF can be derived from [dHb] based on the equations⁷:

[dHb] = (1 - S_{v} O_{2}) \cdot [Hb]

(10)

with [Hb] the haemoglobin concentration in blood, and:

O E F = 1 - \frac{φ \cdot [H b]}{C a O_{2}} \cdot S_{v} O_{2}

(11)

where ϕ is the oxygen binding capacity of haemoglobin (ϕ = 1.34 ml/g) and CaO₂ is the oxygen concentration in arterial blood.

Notably, C and $β_{S E - G E}$ may be modified by vessel topology and intravascular signal contribution to the total BOLD signals even for fixed TEs and readout schemes. These factors necessitate further investigation of the proposed method beyond analytical modelling. We examined the validity of the modelling and the values of the two constants in equation (9) based on Monte Carlo simulations of the BOLD signal (see the Simulations section for further details).

Simulations

Multi-compartmental Monte Carlo modelling of the BOLD signal

Similar to the work of Uludağ et al.,⁴⁴ we implemented a Monte Carlo three-compartment model of the BOLD signal, which comprises venous CBV (CBV_v), capillary CBV (CBV_cap) and tissue volume (1 − CBV_dHb), where CBV_dHb is the sum of CBV_v and CBV_cap. The total MRI signal can then be expressed as:

\begin{array}{l} S_{t o t} = (1 - C B V_{d H b}) \cdot S_{e x} + (1 - I n_{s u p}) \cdot \\ (C B V_{v} \cdot S_{i n, v} + C B V_{c a p} \cdot S_{i n, c a p}) \end{array}

(12)

were S_ex is the extravascular signal and S_in,v and S_in,cap are the intravascular signals for the venous and capillary compartment, respectively. Importantly, since the proposed approach exploits the extravascular BOLD signals, its performance is affected by the extent of intravascular signal, which can be suppressed using tailored sequences. Thus, the parameter In_sup in equation (12) was introduced to depict the level of intravascular suppression (ranging between 0 and 1). As in previous work by our group,⁷ the relation between CBV_dHb, CBV_v and CBV_cap was assumed to be modulated by a parameter ρ, such that:

{\begin{matrix} C B V_{d H b} = ρ \cdot C B V_{c a p} \\ C B V_{v} = (ρ - 1) \cdot C B V_{c a p} \end{matrix}

(13)

Notably, ρ = CBV_dHb/CBV_cap modulates the capillary and venular contribution to the BOLD signal and, as such, it also modulates the average vessel radius within the voxel.

An expected value of ρ = 2 was assumed⁷ with a certain level of variability allowed.

The GE and SE signals for each compartment were assumed equal to:

{\begin{matrix} S_{G E} = S_{o} \cdot e^{- R_{2}^{*} \cdot T E_{G E}} \\ S_{S E} = S_{o} \cdot e^{- R_{2} \cdot T E_{S E}} \end{matrix}

(14)

with TE_GE and TE_SE = 30 and 85 ms, respectively (matching our in-vivo recordings).

Monte Carlo modelling within the ODIN framework⁴⁵ was used to simulate extravascular R_2,ex* and R_2,ex as a function of CBV_dHb and magnetic susceptibility χ for venules and capillaries.⁴⁴ Radiuses of 5 and 24 μm were chosen to reflect the capillary and the venular scale, respectively. The vasculature was modelled as fully permeable to water, randomly oriented, cylinders of infinite length.

To quantify the [dHb] contribution to the extravascular signal, χ was calculated from vessel oxygen saturation, SO₂ as⁴⁴:

Δ χ = 4 π \cdot H c t \cdot Δ χ_{do} \cdot (| {SO}_{2, off} - {SO}_{2} |)

(15)

where Δχ is the susceptibility difference between the vessel and the surrounding tissue, Δχ_do = 0.264 ppm is the susceptibility of fully deoxygenated blood,⁴⁵ Hct is the haematocrit (computed assuming a proportionality with [Hb], Hct/[Hb] = 3%dl/g).⁴⁶ SO_2,off was assumed to be SO_2,off = 0.95.⁴⁵

Random vessel orientations were accounted for by varying the angle between the main field and the cylinder in the range from 0 to π/2 (in 64 linear steps). A sinθ weighting was applied to the signal values to account for the density of possible vessel angles with respect to the static field. For each simulated radius and vessel orientation, 10⁵ random walks were simulated from randomly seeded starting points. The temporal step size of each simulation was 1 ms, with the displacement at each step generated by a gaussian distribution to represent isotropic diffusion.

Diffusion was assumed to be 0.76 μm²/ms, as is typical in GM measurements,^47,48 and a CBV_dHb = 3% was used for explicit simulations and assumed to have a linear effect on relaxation rates when a different blood volume was required. Gradient echo times of 17, 32.5 and 48 ms were used to simulate GE signals which were monoexponentially fitted to calculate R_2,ex* values, while a mono-exponential fit from the initial amplitude and a spin echo at 85 ms was used to estimate R_2,ex.

The total extravascular R₂* and R₂ were calculated from a linear combination of a fixed non-blood relaxation rate arising from the tissue and the extravascular contributions from the venous and capillary blood compartments:

{\begin{matrix} R_{2, e x}^{*} = R_{2, e x, t}^{*} + R_{2, e x, v}^{*} + R_{2, e x, c a p}^{*} \\ R_{2, e x} = R_{2, e x, t} + R_{2, e x, v} + R_{2, e x, c a p} \end{matrix}

(16)

Tissue $R_{2, e x, t}^{*}$ and $R_{2, e x, t}$ were assumed equal to 21 and 13 Hz at 3 T,⁴⁴ respectively. The extravascular GE and SE signals were calculated from R₂*,_ex and R_2,ex according to equation (14).

Intravascular R₂* and R₂ for veins and capillaries were derived from the phenomenological equations reported by Zhao et al. and used to estimate intravascular signals from equation (14).⁴⁹ Capillaries were assumed to have an oxygen saturation midway between arteries and veins. Intravascular signals (in veins and capillaries) were combined with the extravascular signals based on equation (12) to obtain total gradient and spin echo signals (S_{tot_GE} and S_{tot_SE}).

During vasodilation and modulation of CBF, the relative CBV change for the two intravascular compartments was assumed to follow the Grubb relation (α = 0.18),⁴¹ whereas the [dHb] was selectively affected by CBF changes (constant CMRO₂). The BOLD_GE and BOLD_SE signal changes were computed as:

{\begin{matrix} \frac{Δ B O L D_{G E}}{B O L D_{G E}} = \frac{S_{t o t_G E, m} - S_{t o t_G E}}{S_{t o t_G E}} \\ \frac{Δ B O L D_{S E}}{B O L D_{S E}} = \frac{S_{t o t_S E, m} - S_{t o t_S E}}{S_{t o t_S E}} \end{matrix}

(17)

where the subscript tot_GE and tot_SE represent the total gradient and spin echo signals, respectively, and m represents the signal arising from vasodilatory modulation. Notably, from now on, we will refer to the total BOLD_GE and BOLD_SE signal changes (left-hand side of equation (17)) simply as BOLD_GE and BOLD_SE.

Forward and inverse model simulations

To address the method’s sensitivity to confounding factors and the validity of the Davis model, we ran forward and inverse modelling simulations by exploiting the Monte Carlo modelling of the BOLD signal. The physiological parameters affecting the BOLD_GE and BOLD_SE changes, such as [Hb], CaO₂, OEF, flow modulation (ΔCBF/CBF) and ρ were varied as random variables spanning plausible physiological ranges.

The inversion model was implemented using two strategies: (1) using a grid search approach applied to Monte Carlo modelling (OEF: 0–1, 10⁻³ step, fixed ΔCBF/CBF = 20%, CBV_dHb = 3% and ρ = 2) and (2) explicitly inverting the Davis model (equations (9)–(11)), with parameters fitted using forward Monte Carlo simulations. Figure 1 illustrates the main simulated parameters and their corresponding probabilistic distributions. Notably, the varying parameters used to compute the forward model were either assumed to be known (i.e. [Hb], CaO₂, BOLD_GE, BOLD_SE) or were fixed in the inversion model (ΔCBF/CBF, CBV_dHb, ρ), while the remaining were inferred (OEF, SvO₂ and [dHb]).

Figure 1.

Random variables used to perform the multi-compartmental Monte Carlo forward and inverse model simulations. The variables reported in light grey were assumed to be known by the 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.

The effects of intravascular signal suppression (In_sup), variability in ρ and experimental noise on the inversion model performance were also explored. For the latter, white gaussian noise was added to BOLD signal changes, and SNR was calculated as the expected signal change with the increase in blood flow normalised by noise standard deviation.

Participants

Twenty-eight healthy subjects (19 males; age, mean ± standard deviation: 31 ± 8 years) were enrolled, after providing written informed consent. The study was performed in accordance with the Declaration of Helsinki and approved by the Institutional Review Board (protocol number 14-C, 11/12/2023) of the Department of Neurosciences, Imaging and Clinical Sciences, University ‘G d’Annunzio’ of Chieti-Pescara, Italy.

MRI acquisition

Data were acquired at the Institute for Advanced Biomedical Technologies, University ‘G d’Annunzio’ of Chieti-Pescara, Italy, on a 3 T MAGNETOM Prisma scanner (Siemens Healthineers AG, Forchheim, Germany), featuring a 32-channel receive-only head coil. The fMRI measurements were conducted with a customised pCASL research sequence. An BOLD-ASL sequence was developed based on a manufacturer’s product sequence utilising a dual-excitation (DEXI) echo planar imaging (EPI) 2D readout that includes a short TE, TE1 = 10 ms for ASL, and, following a second excitation, a longer TE, TE2 = 30 ms for BOLD_GE.^7,51 The separate readouts enabled optimal background suppression during the ASL image acquisition while allowing a sufficient amount of longitudinal magnetisation recovery for BOLD readout. Furthermore, a refocussing pulse was added after the second excitation for each slice to generate SE weighting, with TE3 = 85 ms. Pre-labelling saturation and background suppression were included.⁵² The ASL labelling duration (τ) and post-labelling delay (PLD) were set to 1.5 s, and GRAPPA acceleration was applied (factor = 3). An effective TR of 5 s was utilised to acquire 14 slices, with in-plane resolution of 3.4 × 3.4 mm² and slice thickness = 7 mm, with a 30% slice gap for full brain coverage.

fMRI data were collected during a BH task. The BH protocol included a baseline period of 60 s followed by 10 cycles of post-expiratory breath-holding, each lasting 20 s, interleaved with a 40 s-normal breathing, totalling 11 min and 28 s.²⁵ Subjects were visually cued with instructions executed using E-Prime 3.0 stimulus presentation software (Psychology Software Tools, Pittsburgh, PA, USA). Participants were directed to maintain a neutral diaphragm position during BH and to completely exhale at the end of each BH.⁵³

During the fMRI recordings, CO₂ partial pressure in the exhaled air was measured using a nasal cannula connected to a gas analyser system (ADInstruments, Dunedin, New Zealand).

Two proton density images (S₀) were obtained for susceptibility distortion correction and ASL calibration with pCASL labelling and background suppression pulses turned off, with TR = 7 s and TE = 10 ms, and opposite phase encoding directions. A MP2RAGE T1-weighted scan was conducted for registration and brain segmentation purposes (matrix 176 × 256 × 256, 1 mm isotropic resolution, TR/TE = 5000/3.58 ms, TI1/TI2 = 700/2500 ms).⁵³ In a subgroup of 16 subjects, blood T₂ in the SSS was estimated from TRUST acquisition (inversion time = 1020 ms, $τ_{C P M G}$ = 10 ms, labelling thickness 100 mm and T₂-prepared effective (e)TEs 0, 40, 80 and 160 ms).¹² Please refer to Table S1 in Supplementary Materials for additional information about MRI sequences.

Data processing

Gas recordings pre-processing

End-tidal O₂ and CO₂ partial pressures (PetO₂ and PetCO₂), taken as a surrogate measure of arterial oxygen (PaO₂) and CO₂ (PaCO₂) partial pressures, were extracted using in-house software in MATLAB (version R2022b).^7,54

The PetO₂ and PetCO₂ traces were calculated by isolating the minima (for the O₂) and the peaks (for CO₂) of the traces from the gas analyser. The end-tidal signals were resampled at the fMRI TR and shifted to account for time lags between expiration and recordings. Baseline PaO₂ and PaCO₂ values were estimated in the first 60 s.

Baseline PaO₂ was used to infer arterial oxygen saturation, SaO₂, through the equation⁴²:

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

(18)

where h is the Hill constant (h = 2.8) and P₅₀ is oxygen pressure when haemoglobin is half saturated (P₄₉ = 26 mmHg). CaO₂ was inferred from⁷:

C a O_{2} = φ \cdot [H b] \cdot S a O_{2} + ε \cdot P a O_{2}

(19)

where ε is the oxygen plasma solubility (ε = 0.0031 ml/mmHg/dl). [Hb] was assumed to be 13 g/dl.

PetCO₂ traces were band-pass filtered (Butterworth digital filter, cut-off times of 10 and 150 s), to be used as regressor to estimate BOLD_GE, BOLD_SE and ASL perfusion signal changes in response to the BH stimulus.

Anatomical MRI and fMRI pre-processing

The MP2RAGE UNI image was employed for tissue segmentation (FAST, FSL) and for warping into MNI space (antsRegistration, SyN, ANTs).^56–58

The pre-processing steps of fMRI data involved data normalisation, motion correction, susceptibility distortion correction and filtering and were performed using FSL, ANTs and in-house MATLAB algorithms.^55,56 The two proton density (S₀) images were corrected for susceptibility distortions (Topup, FSL)⁵⁸ and intensity inhomogeneity (N4biasfieldcorrection, ANTs).⁵⁹ The corrected S₀ image, skull-stripped using FSL BET, was rigidly registered to the T1-weighted image, and the transformation matrix was inverted to bring the GM and WM partial volume estimates into the S₀ space. These were thresholded (th = 0.5) to obtain compartmental masks. An fMRI motion correction pipeline was applied to the fMRI BH data (Supplementary Materials).⁶⁰ The motion corrected volumes were then rigidly registered (ANTs) to the brain-extracted S₀, acquired with the same phase encoding direction as the functional scans. The fMRI volumes for each echo were then corrected for susceptibility distortions (FSL ApplyTopup).

The perfusion signals (∆S) in S₀ space were obtained through surround subtraction of the fMRI timecourses at TE1¹⁰ and converted to CBF in quantitative units of ml/100 g/min through the pCASL single compartment kinetic model¹⁰:

C B F = \frac{6000 \cdot λ \cdot e^{\frac{P L D}{T 1_{b}}}}{2 \cdot η \cdot η_{i n v} \cdot T 1_{b} \cdot (1 - e^{- \frac{τ}{T 1_{b}}})} \cdot (\frac{Δ S}{S_{0}})

(20)

with λ the water partition coefficient (λ = 0.9 ml/g), T1_b the T1 of arterial blood (T1_b = 1.67 s), η the tagging inversion efficiency (η = 0.85) and η_inv a scaling factor accounting for reduction in tagging efficiency due to background suppression (η_inv = 0.88).⁶²

Surround averaging was performed on the BOLD signals to eliminate ASL contamination. CBF and BOLD signals were expressed as relative changes with respect to the baseline. All three fMRI signals were band-pass filtered in analogy with the PetCO₂ trace.

qfBOLD, cfMRI and TRUST data processing

Processing was performed in MATLAB. The evaluation of voxel-wise BOLD and CBF modulation in response to the BH task was performed using linear regression,⁶² where the filtered PetCO₂ trace served as the independent variable. The filtered PetCO₂ was allowed to shift by ±10 s to account for haemodynamic lags.²⁵ The regression provided estimates of cerebrovascular reactivity (CVR, signal change per unit of PetCO₂ change), which were multiplied by a metric of modulation in the PetCO₂ trace (the difference between the 95th percentile and the fifth percentile of the signal) to obtain BOLD and CBF signal changes in response to the BH tasks.

qfBOLD analysis to infer OEF maps was performed based on voxel-wise ratio of the BOLD_SE and the BOLD_GE signal changes in response to the BH task. This data was input into an inversion model based on the Davis model (equations (9)–(11)) which had been fitted to Monte Carlo simulations (please refer to the Results section). For within-session repeatability assessment, the same processing was performed by dividing the filtered traces by considering separately the first five BH and the last five BH.

A cfMRI analysis was also performed based on the comparison of the BOLD_GE and ASL signal modulations in response to BH and by employing a single-calibration framework we developed to infer OEF in the GM.^7,25 The approach is described in detail in Chiarelli et al.,⁷ but essentially infers OEF from the maximum BOLD_GE modulation (M_GE).^7,25 Briefly, using the central volume theorem, the model expresses the baseline CBV_dHb as the product of CBF (that we measure at baseline with ASL) and the mean transit time (MTT) within the compartment of interest to then link the MTT with OEF via the flow-diffusion model of oxygen transport, thus solving the problem of decoupling the dependence of M_GE on CBV_dHb and OEF. The model requires the assumption of a negligible oxygen tension at the mitochondria, which is plausible in the healthy brain. The equation relating M_GE to OEF is^7,25:

\begin{array}{l} M_{G E} = T E_{G E} \cdot \frac{A_{G E} \cdot ρ}{k} \cdot \\ \frac{C B F \cdot O E F \cdot C a O_{2} \cdot {((1 - \frac{C a O_{2}}{φ [H b]} \cdot (1 - O E F)) \cdot [H b])}^{β_{G E}}}{(P_{50} \cdot \sqrt[h]{\frac{2}{O E F} - 1})} \end{array}

(21)

were k is the effective oxygen permeability of brain tissue. Equation (5) was used to estimate M_GE and equation (21) was inverted to estimate the OEF using a grid search approach (OEF between 0 and 1, in steps of 10⁻³). We assigned a value of 1.3 to β_GE and a value of 8.8/s/g^βdl^β/(μmol/mmHg/ml/min) to the term (A_GE ρ)/k, matching our previously established in-vivo estimations.⁷

Global venous T₂ and OEF was also estimated from the TRUST acquisition through the fitting of signal decay with effective TE in the SSS and a calibration model.⁶³

Statistical analysis

Statistical analysis was performed in MATLAB. Root mean square errors (RMSEs) and Pearson’s correlations were evaluated to assess associations between the variables of interest. T-tests were conducted to evaluate statistical significance. A p < 0.05 was considered significant.

Results

Simulations

Figure 2 shows the results of forward model simulations when fixing [Hb], CaO₂ and ρ. Results for vascular suppression of 0%, 50% or 100% are reported. Notably, when [Hb] and CaO₂ are fixed, the OEF is linearly related to [dHb] (from which OEF is derived through equations (10) and (11) and to which BOLD methods are sensitive). Figure 2(a) and (b) show the BOLD_GE and the BOLD_SE signal as a function of OEF. As expected, no significant association between OEF and BOLD signals was found due to variability in baseline CBV_dHb and in ΔCBF/CBF (all other parameters except OEF were fixed). On the contrary, when evaluating the BOLD_SE/BOLD_GE signal changes ratio, a clear monotonic dependence on OEF is visible with reduced effects introduced by baseline CBV_dHb and ΔCBF/CBF. Importantly, the rate of monotonic dependance is higher for stronger intravascular signal suppression. A larger suppression of the intravascular BOLD signal also tends to decrease the BOLD_SE/BOLD_GE signal changes ratio for a given OEF value. Based on average experimental values of the BOLD_SE/BOLD_GE signal changes ratio in GM and average global OEF estimation from TRUST, we inferred that our BOLD_GE/BOLD_SE ASL sequence showed an intravascular suppression of BOLD signals ~50% (In_sup = 0.5 in equation (12)).

Reports forward model simulations as in Figure 2 with an intravascular suppression of 50%, and ρ variability (±3 σ/μ) of 0% (Figure 3(a)) and 60% (Figure 3(b)). As ρ modulates the ratio between capillary and venules, which differently affect BOLD_GE and BOLD_SE, the variability in ρ is reflected in a higher noise introduced in the BOLD_SE/BOLD_GE monotonic relation with OEF. We fitted the Davis model (equation (9)) to the data in Figure 3(a) around expected values of OEF (0.25–0.55) for an [Hb] = 13 g/dl, obtaining C = 0.14 and β_SE-GE = 0.85. Importantly, C and β_SE-GE, among other factors, are influenced by the extent of intravascular suppression.

Figure 2.

Forward model simulations when fixing [Hb], CaO₂ and ρ (ρ = CBV_dHb/CBV_cap). The images show (a) BOLD_GE modulation, (b) BOLD_SE modulation and (c) BOLD_SE/BOLD_GE signal change ratio as a function of OEF, for three levels of intravascular signal suppression (0%, 50% and 100%).

Figure 3.

Forward model simulations when fixing [Hb], CaO₂ and In_sup (In_sup = level of intravascular signal suppression, between 0 and 1). The subplots show BOLD_SE/BOLD_GE signal change ratio as a function of OEF for (a) ρ variability (±3 σ/μ) = 0% and (b) ρ variability = 60%.

Reports the results obtained with Monte Carlo forward modelling and inversion modelling implemented either by inverting the same Monte Carlo modelling or the Davis analytical modelling within the fitting OEF range used. Figure 4(a) shows correlation and Bland–Altman plots obtained with forward model simulations as presented in Figure 1 but with a fixed ρ = 2. Figure 4(b) reports the same simulation, but adding variability to ρ as reported in Figure 1. Intravascular signal suppression was fixed at 50% for both BOLD signals. The inversion models were able to estimate the OEF based on the BOLD_SE/BOLD_GE signal changes ratio with an OEF RMSE of around 2% and 3.5%, respectively. Notably, the low RMSE obtained from inverting the Davis model fitted to the Monte Carlo simulations corroborates the validity of the analytical framework presented in the “Methods” section (equation (9)).

Figure 4.

Correlation (left image) and Bland–Altman plot (right image) of the estimated OEF versus the simulated OEF using Monte Carlo modelling (black dots) or Davis analytical modelling (blue dots) for inversion and Monte Carlo forward modelling with In_sup = 50% (In_sup = level of intravascular signal suppression, between 0 and 1) and (a) fixing ρ = 2 (ρ = CBV_dHb/CBV_cap) or (b) considering ρ as a random variable with a variability of 60% (±3 σ/μ). OEF RMSE of Davis analytical modelling as a function of (c) In_sup and ρ variability, and as a function of (d) BOLD_GE SNR and BOLD_SE SNR, presented both as a bi-dimensional image (left image) and as a plot (right image). The arrows in the right images of subplots (c) and (d) depict increasing value of In_sup and BOLD_SE SNR, respectively.

Figure 4(c) reports the OEF RMSE of the inversion using the Davis model as a function of In_sup and ρ variability. OEF RMSE increases as In_sup decreases and ρ variability increases, with the greater effect of ρ variability. With ρ variability up to 100%, the method delivers a RMSE <5%. Figure 4(d) reports the OEF RMSE as a function of BOLD_GE and BOLD_SE SNR for In_sup = 50% and ρ variability = 60%. The RMSE due to signal noise is below 5% when the SNR is above 8 for both BOLD_GE and BOLD_SE modulations.

In vivo data

The average PetO₂ and PetCO₂ at rest were 112.7 ± 4.8 and 34.0 ± 3.0 mmHg, respectively. The BH task was successfully performed by all participants and induced a consistent modulation in PetCO₂, reflecting in a modulation of CBF and BOLD signals.

Figure 5(a) shows the average unfiltered (left) and filtered (right) PetCO₂ traces (top row). The PetCO₂ modulation (difference between the signal 95th and the fifth percentiles) was 6.9 ± 1.9 mmHg. Figure 5(a) also shows the average modulations in the GM for CBF and for BOLD_SE and BOLD_GE, respectively (bottom row). The average GM modulation in CBF was 26% ± 16%, whereas the modulations in the BOLD signals were BOLD_GE = 1.50% ± 0.34% and BOLD_SE = 0.79% ± 0.24%. Figure 5(b) shows regional maps of the BOLD_GE and BOLD_SE modulations (top row), BOLD_GE and BOLD_SE SNR (middle row) and BOLD_SE/BOLD_GE signal change ratio (bottom row, left image), for an exemplar study participant. The distribution of BOLD_SE/BOLD_GE signal changes ratio in GM and WM is also reported (bottom row, right image). Figure 5(c) reports the same maps/plots as in Figure 5(b) but as average values across participants (with across subjects’ t-score maps substituting SNR). The subjects’ average BOLD_SE/BOLD_GE signal changes ratio was 58.8% ± 21.9% in the GM and 74.5% ± 25.5% in the WM.

Figure 5.

(a) Average raw PetCO₂ trace across participants and filtered PetCO₂ (top row) as well as average and filtered CBF and BOLD_GE and BOLD_SE modulations in the GM (bottom row), (b) representative maps for a participant of the study and (c) average maps in MNI space of BOLD_GE and BOLD_SE modulations (top row), of BOLD_GE and BOLD_SE SNR or t-score (middle row) and of BOLD_SE/BOLD_GE signal changes ratio (bottom row, left image). The BOLD_SE/BOLD_GE signal changes ratio distributions in GM and WM are also reported (bottom row, right image).

Figure 6(a) shows a representative OEF map for a participant of the study computed from the BOLD_SE/BOLD_GE signal changes ratio. CBF maps (computed with ASL at baseline) and CMRO₂ maps are also displayed. Histograms of OEF values obtained in the GM and WM are also reported. Figure 6(b) reports the same maps/plots as in Figure 6(a) but as average values across participants in MNI space. Repeatability analysis of OEF maps delivered good voxelwise spatial correlations in GM and WM (across subjects’ r = 0.41 ± 0.09 for GM and r = 0.38 ± 0.09 for WM, respectively, both p’s < 10⁻³) and good-to-excellent correlations on a global basis (r = 0.84 for GM and r = 0.72 for WM, respectively, both p’s < 10⁻³, refer to Figure S1 in Supplementary Materials).

Figure 6.

(a) Representative images for a participant of the study and (b) average maps in MNI space for OEF (distribution values in GM and WM are also reported), CBF and CMRO₂.

Figure 7 reports scatterplots and Bland–Altman plots comparing global OEF values in the GM and WM obtained with qfBOLD with alternative approaches. The average OEF for qfBOLD was OEF = 37.0% ± 2.9% in the GM and OEF = 41.6% ± 2.9% in the WM with a significant correlation (r = 0.88, p < 10⁻³) and different average values between the two compartments (p < 10⁻³). Figure 6(a) compares the proposed method with the single calibration fMRI approach in the GM. The average GM OEF was OEF = 36.6% ± 2.0% for single calibration fMRI. A correlation of r = 0.71 (p < 10⁻³) was obtained in the GM with no bias between the approaches. Figure 7(b) and (c) compare the GM and WM OEF extracted with qfBOLD with the OEF evaluated in the SSS with TRUST in a subset of 16 subjects. The average OEF from TRUST was 40.1% ± 5.2%. The correlation coefficients in the GM and WM were r = 0.51 (p < 0.05) and r = 0.61 (p < 0.01), respectively. TRUST showed slightly larger OEF in the GM and smaller OEF in the WM compared to those obtained using the qfBOLD method (p < 0.05).

Figure 7.

Scatterplots and Bland–Altmann plots comparing extracted global OEF values in the GM and WM with alternative approaches: (a) qfBOLD versus single calibration fMRI for GM and (b) TRUST versus qfBOLD in GM and (c) TRUST versus qfBOLD in WM.

Discussion

We have introduced a novel framework for mapping OEF at rest that relies on comparing (taking the ratio of) functional modulations in the gradient-echo and spin-echo BOLD signals (BOLD_GE and BOLD_SE), both of which are concurrently acquired during isometabolic vasodilation. This method, termed quantitative functional BOLD (qfBOLD), differs from traditional qBOLD,¹⁵ in both its approach and acquisition strategy. Traditional qBOLD quantifies OEF by comparing multi-TE GE and SE acquisitions, often relying on absolute estimations of T₂*, T₂ and T₂′. In contrast, qfBOLD derives OEF from temporal fluctuations in T₂* and T₂, each estimated from BOLD_GE and BOLD_SE acquired at a single TE. As dHb is the only magnetic substance whose concentration oscillates over time, qfBOLD retains the advantage of selective sensitivity to baseline dHb as in cfMRI while offering several improvements.21,24 qfBOLD can decouple OEF from CBV_dHb with a single modulation of brain physiology and without additional assumptions, while eliminating the need for concurrent functional ASL recordings. This absence of ASL potentially greatly enhances SNR and may enable significant improvements in both temporal and spatial resolution, while potentially allowing extension of the method to WM.

Monte Carlo simulations of the BOLD_GE and BOLD_SE signal changes generated within a voxel with randomly oriented microvasculature during increases in CBF demonstrate the feasibility of the qfBOLD approach with the possibility to use the simplistic Davis model within the framework (Figures 1 –4). Notably, the simulated BOLD_SE/BOLD_GE signal changes ratio reveals a monotonic dependence on OEF (Figure 2(c)), which actually arises from a monotonic dependence on dHb levels in blood ([dHb], with a fixed concentration of haemoglobin in blood, [Hb] and arterial oxygen content, CaO₂). Crucially, this relationship is nearly independent of baseline CBV_dHb and modulations in CBF (Figures 2(c) and 3), supporting the key features of the method, that is, its ability to decouple [dHb] and CBV_dHb relying only on vasodilation and without requiring functional CBF measures. The slope of this monotonic relationship is particularly pronounced under conditions of greater intravascular signal suppression.

Modelling and in vivo studies suggest that the intravascular signal accounts for around 30% of the BOLD_GE signal and 50% of the BOLD_SE signal,^27-31 modifying the total signal dependence on [dHb] compared to pure extravascular effects,²⁶ which primarily generate the method’s sensitivity to OEF. Clearly, the intravascular contribution can be modelled and accounted for, or alternatively, it can be suppressed to enhance the method’s sensitivity to OEF, with the added advantage of partially eliminating the effects of large vessels that the method ignores.

Indeed, the monotonic dependence of the signal ratio is also influenced by microvascular topology, as illustrated in Figure 3, where variability in the parameter ρ adds noise to the monotonic relationship between OEF and the BOLD_SE/BOLD_GE ratio (Figure 3(b)). Nonetheless, the method’s ability to infer OEF is demonstrated in Figure 4, where inversion models achieved an RMSE in OEF estimation below 5% (Figure 4(a) and (b)). The inversion of the Davis model delivered sufficiently small RMSE up to ρ variability of around 100% (RMSE always below 5%) for almost any level of intravascular suppression (Figure 4(c)). This performance was contingent on maintaining a SNR >8 for both BOLD_GE and BOLD_SE signal changes following CBF modulation (Figure 4(d)), which is achievable in vivo (see Figure 5(b)). With respect to in vivo applications, the fitted Davis model was preferred for modelling inversion compared to the full Monte Carlo model, as simpler models are preferable in the presence of experimental noise.

The qfBOLD method exploiting the Davis model inversion was compared in vivo to a single cfMRI approach that we recently developed based on BOLD_GE/ASL recordings acquired during BH.^7,25 We acquired fMRI data during a BH task using a sequence that integrates functional BOLD_GE, BOLD_SE and ASL. Additionally, the new method was also compared to global OEF measures in the sagittal sinus using a relaxometry-based method (TRUST).^12,63

By comparing the average BOLD_SE/BOLD_GE signal changes ratio in GM during the hypercapnic task (Figure 5(c)) with simulations (Figures 2 –4), we determined that the developed sequence achieved ~50% intravascular suppression of the BOLD signals, consistent with an average OEF of ~40%. This estimated suppression is partly due to T1 weighting resulting from the short time interval between ASL and BOLD readouts, combined with the effects of gradients introduced around the spin-echo refocussing pulse. Moreover, the intravascular suppression model parameter may account for inaccuracies in the intravascular signal modelling, which relies on phenomenological equations.⁴⁹

While the OEF maps generated by inverting the developed model appear plausible (Figure 6), it is important to note that macrovascular contamination is present in the maps, with clear regions where the BOLD_SE signal tends to be smaller than the BOLD_GE signal, leading to an underestimation of OEF. To mitigate the influence of errors induced by large vessels we focussed exclusively on voxels with a BOLD_SE/BOLD_GE signal changes ratio above 30% when computing global values for correlation analysis with alternative methods.

The qfBOLD method yielded an average OEF = 37.0% ± 2.9% in GM, compatible with the single calibration fMRI approach (OEF = 36.6% ± 2.0%) and the TRUST method (OEF = 40.1% ± 5.2%). OEF values showed good within-session repeatability at a voxelwise level (spatial correlations: r = 0.41 ± 0.09 in GM and r = 0.38 ± 0.09 in WM, both p’s < 10⁻³), and good-to-excellent repeatability at a global level (r = 0.84 in GM and r = 0.72 in WM, p’s < 10⁻³). Importantly, OEF values measured with qfBOLD showed significant correlations with those estimated with the former alternative approaches on a global basis for GM (qfBOLD vs cfMRI r = 0.71, p < 10⁻³, qfBOLD vs TRUST r = 0.51, p < 0.05) and WM (qfBOLD vs TRUST r = 0.61, p < 0.01). WM exhibited higher average OEF (OEF = 41.6% ± 2.9%) compared to GM, probably induced by modelling limitations (see below). The global correlations we obtained between the OEF estimated with qfBOLD and that obtained through alternative approaches such as cfMRI and TRUST, strongly suggest the validity of the newly proposed approach.

However, several important assumptions in the vascular modelling warrant consideration. The assumption of isotropic vessel orientation, while predominantly valid in GM, may not hold in WM, where complex topology and anisotropic vasculature could contribute to an overestimation of OEF.^64–66 Additionally, restricted diffusion in the extravascular space of WM may affect measurement accuracy, as water diffusion was assumed fixed in the model and tuned to GM. Finally, although the method focusses on modulations in T₂ and T₂* due to dHb, making it less susceptible than qBOLD to additional sources of magnetic susceptibility, it may retain residual sensitivity to WM tissue anisotropy and susceptibility-induced signal decay modifications, for example, those related to myelin.

While our results correlate well with two alternative methods further development is necessary. Future advancements in the method would require additional modelling and data acquisition improvements. Modelling may be refined by integrating realistic topologies of vascular compartments and variability in water diffusion, particularly for exploring WM. Experimentally, by eliminating ASL, the focus should be on acquiring high-quality BOLD_GE and BOLD_SE data with high spatiotemporal resolution and on maximising SNR for optimal method stability and repeatability. Moreover, macrovascular signal suppression and intravascular suppression of microvascular signals are warranted to reduce bias and variance.

The method repeatability between sessions and possibly MRI scanners should be assessed and, importantly, the approach should be tested in diseases where larger variability in vascular topology may decrease sensitivity to OEF. In addition, similarly to cfMRI approaches, the method requires a vascular reserve, which may be absent in diseases with compromised vasculature. In practice, when hypercapnia is induced via breath-holding, the task requires subject compliance.²⁵ Another practical limitation of the approach, particularly for clinical settings, pertains to the need for acquiring end-tidal partial pressures. We used PetCO₂ measures to derive the hypercapnic modulatory signal, and the PetO₂ signal to derive baseline parameters of arterial oxygenation. However, an alternative approach that avoids gas recordings was previously presented by our group for cfMRI and it can be applied to qfBOLD. The approach is based on deriving the vasodilatory signal from the global BOLD signal.²⁵ Moreover, baseline oxygenation may be fixed to standard values, as the error introduced in the OEF estimates is within a few percentage points.⁴² Nonetheless, hypercapnia itself may also introduce a bias in the measurements, as the assumption of isometabolism may not completely hold; some studies report that hypercapnia mildly reduces neural activity and CMRO₂.^42,65–67 However, this effect should introduce a smaller bias in the qfBOLD method compared to cfMRI, as changes in CMRO₂ would similarly affect both BOLD_SE and BOLD_GE.

The new qfBOLD method appears viable and holds promise for accurate mapping of brain oxygen extraction and consumption with MRI, offering advantages over current alternative methods such as qBOLD, QSM and cfMRI.

Conclusion

We introduced a novel framework for OEF mapping with MRI, termed quantitative functional BOLD (qfBOLD). qfBOLD leverages the extravascular temporal dynamics of gradient-echo BOLD and spin-echo BOLD signals during isometabolic vasodilation. This innovative approach enhances sensitivity to baseline dHb compared to qBOLD. Furthermore, it effectively decouples OEF from CBV_dHb with a single modulation in brain physiology, without the need for concurrent ASL recordings, which are necessary in calibrated fMRI approaches. The avoidance of ASL improves SNR and allows for increased spatiotemporal resolution, enabling its application also in WM. While promising, the method’s reliance on specific vascular modelling assumptions, particularly in WM, highlights the need for further refinement in future applications. Ultimately, the qfBOLD method shows substantial potential for providing accurate and reliable assessments of brain oxidative metabolism. Further studies will be essential to optimise the method, for example, by maximising intravascular signal suppression, and explore its applicability across various pathological conditions.

Supplemental Material

sj-docx-1-jcb-10.1177_0271678X261453807 – Supplemental material for Quantitative functional BOLD (qfBOLD): A combined gradient-echo and spin-echo framework for oxygen extraction fraction (OEF) mapping with functional MRI

Supplemental material, sj-docx-1-jcb-10.1177_0271678X261453807 for Quantitative functional BOLD (qfBOLD): A combined gradient-echo and spin-echo framework for oxygen extraction fraction (OEF) mapping with functional MRI by Antonio Maria Chiarelli, Lucie Chalet, Sara Pomante, Davide Di Censo, Alessandra Caporale, Emma Biondetti, Fabrizio Fasano, Domenico Zaca, Giulia Rocco, Manuela Carriero, Francesca Graziano, Elizabeth Jane Fear, Maria Eugenia Caligiuri, Richard Geoffrey Wise and Michael Germuska in Journal of Cerebral Blood Flow & Metabolism

Footnotes

Author contributions

AMC and MG developed imaging and analysis methods and analysed and interpreted the data. SP, DDC, LC and AC analysed the data. FF and DZ helped develop the MRI sequence. AMC and MG drafted the manuscript. AMC, MG, MEC and RGW conceived the project. EJF, FG, MC, LC and GR set up and executed the experiment. LC, AC and EB contributed to data analysis and manuscript revision. All authors revised and approved the final submission.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: EU, Italian MUR, PNR and PRIN, n. 2022BERM2F, n. 2022MHMSSJ. EU, Italian MUR, PNRR, n. P20225AEEE, n. P2022ESHT4, n. ECS00000041-VITALITY, n. PE0000006-MNESYS. EB has received funding from the EU under the Marie Skłodowska-Curie Grant Agreement n. 101066055.

MG is supported by the Wellcome Trust (220575/Z/20/Z) and received funding from the Engineering and Physical Sciences Research Council (EP/S025901/1).

Declaration of conflicting interests

The authors declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: Fabrizio Fasano and Domenico Zaca work for Siemens Healthineers, the producer of the scanner used to acquire the data presented in the study. The other authors declare no competing financial and non-financial interests.

ORCID iDs

Antonio Maria Chiarelli

Lucie Chalet

Emma Biondetti

Data availability statement

Data and coding is available upon reasonable request.

Supplemental material

Supplemental material for this article is available online.

References

Senior

AE.

ATP synthesis by oxidative phosphorylation. Physiol Rev 1988; 68: 177–231.

Magistretti

Allaman

Brain energy metabolism. In: Pfaff

(ed.) Neuroscience in the 21st century. New York, NY: Springer, 2013, pp.1591–1620.

Lassen

Christensen

Physiology of cerebral blood flow. Br J Anaesth 1976; 48: 719–734.

Gjedde

The pathway for oxygen in brain. APMIS Suppl 2003; 2003: 146–153.

Gjedde

Poulsen

Østergaard

. On the oxygenation of hemoglobin in the human brain. In: Eke

Delpy

(eds.) Oxygen transport to tissue XXI. Boston, MA: Springer, 1999, pp.67–81.

Habler

Messmer

KFW

. The physiology of oxygen transport. Transfus Sci 1997; 18: 425–435.

Chiarelli

Germuska

Chandler

, et al. A flow-diffusion model of oxygen transport for quantitative mapping of cerebral metabolic rate of oxygen (CMRO₂) with single gas calibrated fMRI. J Cereb Blood Flow Metab 2022; 42: 1192–1209.

Rodgers

Detre

Wehrli

FW.

MRI-based methods for quantification of the cerebral metabolic rate of oxygen. J Cereb Blood Flow Metab 2016; 36: 1165–1185.

Fan

Moradi

, et al. Quantification of brain oxygen extraction and metabolism with [¹⁵O]-gas PET: a technical review in the era of PET/MRI. Neuroimage 2020; 220: 117136.

10.

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.

11.

Biondetti

Cho

Lee

Cerebral oxygen metabolism from MRI susceptibility. Neuroimage 2023; 276: 120189.

12.

Quantitative evaluation of oxygenation in venous vessels using T₂-relaxation-under-spin-tagging MRI. Magn Reson Med 2008; 60: 357–363.

13.

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.

14.

Liu

Haist

Dubowitz

, et al. Cerebral blood volume changes during the BOLD post-stimulus undershoot measured with a combined normoxia/hyperoxia method. Neuroimage 2019; 185: 154–163.

15.

Zhu

Yablonskiy

DA.

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

16.

Zhang

Liu

Gupta

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

17.

Hametner

Endmayr

Deistung

, et al. The influence of brain iron and myelin on magnetic susceptibility and effective transverse relaxation – a biochemical and histological validation study. Neuroimage 2018; 179: 117–133.

18.

Duyn

Schenck

Contributions to magnetic susceptibility of brain tissue. NMR Biomed 2017; 30: e3546.

19.

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.

20.

Hoge

RD.

Calibrated fMRI. Neuroimage 2012; 62: 930–937.

21.

Germuska

Wise

RG.

Calibrated fMRI for mapping absolute CMRO₂: practicalities and prospects. Neuroimage 2019; 187: 145–153.

22.

Bulte

Kelly

Germuska

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

23.

Germuska

Bulte

DP.

MRI measurement of oxygen extraction fraction, mean vessel size and cerebral blood volume using serial hyperoxia and hypercapnia. Neuroimage 2014; 92: 132–142.

24.

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.

25.

Driver

Chiarelli

Chandler

, et al. Breath-hold calibrated fMRI mapping of absolute cerebral metabolic rate of oxygen metabolism (CMRO₂): an assessment of the accuracy and repeatability in a healthy adult population. Imaging Neurosci 2024; 2: 1–14.

26.

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.

27.

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.

28.

Boxerman

Hamberg

Rosen

, et al. MR contrast due to intravascular magnetic susceptibility perturbations. Magn Reson Med 1995; 34: 555–566.

29.

Fujita

Extravascular contribution of blood oxygenation level-dependent signal changes: a numerical analysis based on a vascular network model. Magn Reson Med 2001; 46: 723–734.

30.

Jochimsen

Norris

Mildner

, et al. Quantifying the intra- and extravascular contributions to spin-echo fMRI at 3 T. Magn Reson Med 2004; 52: 724–732.

31.

Zhong

Kennan

Fulbright

, et al. Quantification of intravascular and extravascular contributions to BOLD effects induced by alteration in oxygenation or intravascular contrast agents. Magn Reson Med 1998; 40: 526–536.

32.

Buxton

RB.

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

33.

Liu

van Zijl

PCM

, et al. Three-dimensional whole-brain mapping of cerebral blood volume and venous cerebral blood volume using Fourier transform-based velocity-selective pulse trains. Magn Reson Med 2021; 86: 1420–1433.

34.

Gould

Tsai

Kleinfeld

, et al. The capillary bed offers the largest hemodynamic resistance to the cortical blood supply. J Cereb Blood Flow Metab 2017; 37: 52–68.

35.

Brittain

McCabe

Khatun

, et al. An MRI-histological study of white matter in stroke-free SHRSP. J Cereb Blood Flow Metab 2013; 33: 760–763.

36.

Troprès

Grimault

Vaeth

, et al. Vessel size imaging. Magn Reson Med 2001; 45: 397–408.

37.

Shen

Ahearn

Clemence

, et al. Magnetic resonance imaging of the mean venous vessel size in the human brain using transient hyperoxia. Neuroimage 2011; 55: 1063–1067.

38.

Lemasson

Valable

Farion

, et al. In vivo imaging of vessel diameter, size, and density: a comparative study between MRI and histology. Magn Reson Med 2013; 69: 18–26.

39.

Germuska

Meakin

Bulte

DP.

The influence of noise on bold-mediated vessel size imaging analysis methods. J Cereb Blood Flow Metab 2013; 33: 1857–1863.

40.

Grubb

Raichle

Eichling

, et al. The effects of changes in PaCO₂ cerebral blood volume, blood flow, and vascular mean transit time. Stroke 1974; 5: 630–639.

41.

Leung

Tachtsidis

Tisdall

, et al. Estimating a modified Grubb’s exponent in healthy human brains with near infrared spectroscopy and transcranial Doppler. Physiol Meas 2009; 30: 1–12.

42.

Chiarelli

Germuska

Di Censo

, et al. Multiparametric mapping of brain oxygen consumption with resting state calibrated functional MRI. Neuroimage 2025; 320: 121465.

43.

Chiacchiaretta

Romani

Ferretti

Sensitivity of BOLD response to increasing visual contrast: spin echo versus gradient echo EPI. Neuroimage 2013; 82: 35–43.

44.

Uludağ

Müller-Bierl

Uğurbil

An integrative model for neuronal activity-induced signal changes for gradient and spin echo functional imaging. Neuroimage 2009; 48: 150–165.

45.

Jochimsen

von Mengershausen

ODIN – object-oriented development interface for NMR. J Magn Reson 2004; 170: 67–78.

46.

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.

47.

Jensen

Chandra

Strong field behavior of the NMR signal from magnetically heterogeneous tissues. Magn Reson Med 2000; 43: 226–236.

48.

Bihan

Warach

SJ.

Diffusion and perfusion magnetic resonance imaging: applications to functional MRI. J Comput Assist Tomogr 1995; 19: 844.

49.

Zhao

Clingman

Närväinen

, et al. Oxygenation and hematocrit dependence of transverse relaxation rates of blood at 3T. Magn Reson Med 2007; 58: 592–597.

50.

Schmithorst

Hernandez-Garcia

Vannest

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

51.

Okell

Chappell

Kelly

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

52.

Murphy

Harris

Wise

RG.

Robustly measuring vascular reactivity differences with breath-hold: normalising stimulus-evoked and resting state BOLD fMRI data. Neuroimage 2011; 54: 369–379.

53.

Marques

Kober

Krueger

, et al. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. Neuroimage 2010; 49: 1271–1281.

54.

Takano

Sakamoto

Kiyofuji

, et al. A comparison of the end-tidal CO₂ measured by portable capnometer and the arterial PCO₂ in spontaneously breathing patients. Respir Med 2003; 97: 476–481.

55.

Tustison

Cook

Holbrook

, et al. The ANTsX ecosystem for quantitative biological and medical imaging. Sci Rep 2021; 11: 9068.

56.

Jenkinson

Beckmann

Behrens

TEJ

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

57.

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.

58.

Andersson

JLR

Skare

Ashburner

How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. Neuroimage 2003; 20: 870–888.

59.

Tustison

Avants

Cook

, et al. N4ITK: improved N3 bias correction. IEEE Trans Med Imaging 2010; 29: 1310–1320.

60.

Wang

Aguirre

Rao

, et al. Empirical optimization of ASL data analysis using an ASL data processing toolbox: ASLtbx. Magn Reson Imaging 2008; 26: 261–269.

61.

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.

62.

Penny

Friston

Ashburner

, et al. Statistical parametric mapping: the analysis of functional brain images. London, UK: Elsevier, 2011.

63.

Grgac

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

64.

Nonaka

Akima

Hatori

, et al. Microvasculature of the human cerebral white matter: Arteries of the deep white matter. Neuropathology 2003; 23: 111–118.

65.

Schilling

Newton

Tax

CMW

, et al. The relationship of white matter tract orientation to vascular geometry in the human brain. Sci Rep 2025; 15: 18396.

66.

Doucette

Wei

Hernández-Torres

, et al. Rapid solution of the Bloch-Torrey equation in anisotropic tissue: application to dynamic susceptibility contrast MRI of cerebral white matter. NeuroImage 2019; 185: 198–207.

67.

Baas

KPA

Shen

, et al. Venous blood oxygenation measurements using TRUST and T₂-TRIR MRI during hypoxic and hypercapnic gas challenges. J Magn Reson Imaging 2023; 58: 1903–1914.

68.

Deckers

Bhogal

Dijsselhof

, et al. Hemodynamic and metabolic changes during hypercapnia with normoxia and hyperoxia using pCASL and TRUST MRI in healthy adults. J Cereb Blood Flow Metab 2022; 42: 861–875.

69.

James

Sanggaard

Akif

, et al. Spatiotemporal features of neurovascular (un)coupling with stimulus-induced activity and hypercapnia challenge in cerebral cortex and olfactory bulb. J Cereb Blood Flow Metab 2023; 43: 1891–1904.

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