## Abstract

Magnetic resonance imaging (MRI) offers the possibility to non-invasively map the rate of cerebral metabolic oxygen consumption (CMRO_{2}), which is essential for understanding and monitoring neural function in both health and disease. Existing methods of mapping CMRO_{2}, based on respiratory modulation of arterial spin labelling (ASL) and blood oxygen level dependent (BOLD) signals, require lengthy acquisitions and independent modulation of both arterial oxygen and carbon dioxide levels. Here, we present a new simplified method for mapping the rate of cerebral oxygen metabolism that can be performed using a simple breath-holding paradigm. The method incorporates flow-diffusion modelling of oxygen transport and physiological constraints to create a non-linear mapping between the maximum BOLD signal, M, baseline blood flow (CBF_{0}), and CMRO_{2}. A gradient boosted decision tree is used to learn this mapping directly from simulated MRI data. Modelling studies demonstrate that the proposed method is robust to variation in cerebral physiology and metabolism. This new gas-free methodology offers a rapid and pragmatic alternative to existing dual-calibrated methods, removing the need for specialist respiratory equipment and long acquisition times. In-vivo testing of the method, using an 8-minute 45 second protocol of repeated breath-holding, was performed on 15 healthy volunteers, producing quantitative maps of cerebral blood flow (CBF), oxygen extraction fraction (OEF), and CMRO_{2}.

## 1 Introduction

A continuous supply of oxygen to the brain is essential for life and any restriction of this supply can have significant consequences for individuals ^{1}. The current gold standard for mapping the cerebral metabolic rate of oxygen consumption (CMRO_{2}) is 15-oxygen labelled positron emission tomography ^{2}. However, this method has substantial drawbacks, including long acquisition times, the use of ionizing radiation, and the need for local production of 15-oxygen labelled tracers. Due to these limitations, there is great interest in developing rapid and non-invasive methods that can be used to safely and quickly map CMRO_{2}. In recent years, a number of MRI methods ^{3–6} has been proposed as an alternative to PET based methods. One promising approach is the so-called dual-calibrated fMRI method ^{7–9}. This approach combines biophysical modelling of the blood oxygen level dependent (BOLD) signal, and modulation of cerebral blood flow and blood oxygenation with controlled respiratory stimuli. This approach has been successfully applied in disease ^{10, 11}, sports related head impact ^{12}, and pharmacological modulation studies ^{13}. However, the dual-calibrated method requires specialist respiratory equipment and experienced operators to acquire the data. These requirements limit the wider adoption of such a method both in research and clinical settings. In this work we present a straightforward gas-free approach for estimating the resting cerebral metabolic rate of oxygen consumption (CMRO_{2,0}). The method uses a simple repeated breath-holding paradigm to provide robust estimates of CMRO_{2,0}. The data can be acquired in less than 10 minutes with minimal operator training.

The dual-calibrated method can be conceptualised as making two independent measurements of M (the maximum possible BOLD signal); one measurement is dependent on the oxygen extraction fraction (OEF), and one is not ^{7}. OEF is then estimated by assigning a value that brings the two measurements of M into agreement. In our proposed method only one measurement of M is required. Estimates of M can be acquired with a calibrated-fMRI experiment using either a hypercapnic gas challenge ^{14} or a breath-holding protocol ^{15}, or by transforming a measurement of R_{2}’ ^{16}. In the work presented here we have implemented the modelling and analysis for a repeated breath-holding protocol, allowing mapping of CMRO_{2,0} from a simple imaging protocol without the need for administration of modified gas mixtures.

The proposed method is based on a new formulation of M, which includes a simple model of oxygen exchange ^{17–19}, and substitutes the deoxyhaemoglobin sensitive blood volume (CBV_{v}) for an appropriately scaled capillary blood volume (CBV_{cap}). By combining the newly formulated expression of M with reasonable physiological constraints, we can create a non-linear mapping from M and the resting cerebral blood flow (CBF_{0}) to CMRO_{2,0}. For in-vivo data analysis we employ a machine learning methodology that estimates CMRO_{2,0} directly from the MRI data, without estimating the intermediate parameters. We have previously shown that such an approach, trained on simulated data, is able to make robust estimates of physiological and metabolic parameters when applied to the dual-calibrated fMRI method ^{19}.

We adopted the same concept as used in our recent publication ^{19} and estimate CMRO_{2,0} from a Fourier transformed representation of the MRI time series data. In the proposed implementation we use a popular gradient boosted decision tree algorithm, LightGBM^{20}, to learn the mapping between the simulated MRI data and CMRO_{2,0}. The method is tested via further simulations and in-vivo data acquired from with 15 healthy volunteers.

## 2 MRI Data Modelling

Modelling of the steady-state BOLD signal shows that the rate of transverse relaxation due to deoxyhaemoglobin can be expressed as shown in equation 1 ^{14, 21}.

Where is the transverse relaxation rate due to deoxyhaemoglobin. CBV_{v} is the fractional BOLD sensitive blood volume, which is composed of the venous blood volumes, and capillary blood volumes that are exchanging oxygen with the tissue. β is a field strength dependent constant, and A is a proportionality constant related to the field strength and vessel geometry (in units of s^{−1}g-^{β}dL^{β}). [Hb] is the haemoglobin concentration (g/dL) and SvO_{2} is the draining venous oxygen saturation.

The maximum BOLD signal, M, is obtained simply by multiplying by the acquisition echo time, TE ^{14}.

Using this definition of M, the generic BOLD signal model that describes the signal change due to modulation of arterial oxygen content (CaO_{2}), CBF, or CMRO_{2} can be expressed as shown in equation 3a ^{19}.

Where, Δ*BOLD/BOLD*_{0} is the fractional change in BOLD signal, φ is the oxygen binding capacity of haemoglobin (1.34 ml/g), and the subscript 0 denotes the resting value.

In agreement with the previous literature, we assume that brief periods of breath-holding are iso-metabolic ^{15, 22, 23}. Additionally, we follow recent QSM modelling ^{24} and vascular measurements ^{25} in setting the Grubb exponent, α, to zero. Thus, for a breath-holding stimulus, the generalised BOLD model of equation 3a is simplified to equation 3b. Unlike the standard hypercapnic calibration equation ^{14}, we include the arterial oxygen content (CaO_{2}) as this is known to change during breath-holding ^{26}. Because of this, the equation also maintains a dependence on CMRO_{2,0}.

Typically, the influence of CaO_{2} on the BOLD signal is ignored in calibration experiments based on breath-holding ^{15, 22, 27}. However, here it is included for a more complete description of the BOLD signal. Upon completion of a breath-holding experiment equation 3b will have two unknowns, M and CMRO_{2,0}. Thus, it initially appears that we cannot use a breath-holding paradigm to solve for either parameter. A similar problem is solved in the dual-calibrated methodology by including a hyperoxic measurement, which is predominately sensitive to CBV_{v} ^{28}. Here we resolve the problem by employing oxygen diffusion modelling and physiological constraints during fitting.

Following the modelling of ^{17, 18, 29} we can express the resting cerebral metabolic rate of oxygen metabolism as shown in equation 4a.

Where CBV_{cap} is the capillary blood volume that exchanges oxygen with the tissue, κ is the effective permeability of capillary endothelium and brain tissue (μmol/mmHg/ml/min), P50 is the blood oxygen tension at which haemoglobin is 50% saturated, h is the Hill coefficient (equal to 2.8) and PminO_{2} is the minimum oxygen tension at the mitochondria. In the modelling we calculate the value of P50 from a measure of end-tidal carbon dioxide tension and we assume a fixed value for κ of 3 μmol/mmHg/ml/min ^{19}.

The capillary blood volume exchanging oxygen with the tissue, CBV_{cap}, is a fraction of the BOLD sensitive blood volume (CBV_{v}), i.e. *CBV*_{ν} = *ρ* · *CBV*_{cap}. Thus, by re-arranging equation 4a in terms of CBV_{v} (equation 4b) and substituting into the definition of M we can derive an alternative definition of M as expressed by equation 5.

This definition of M effectively replaces the unknown parameter CBV_{v} with two measurable parameters, CaO_{2,0} and CBF_{0} and two unknown parameters, ρ and PminO_{2}. This can be seen by examining equation 5 and expanding CMRO_{2,0} via the Fick principle, *CMRO*_{2,0} = *CaO*_{2,0} · *CBF*_{0} · *OEF*_{0}. The standard definition of M contains TE, A, [Hb], and SvO_{2}, thus there is already an implicit dependence on OEF_{0} (assuming arterial blood is close to, or fully saturated). The Hill coefficient and the effective permeability of brain tissue can be assumed constant, and the value of P50 can be estimated from the end-tidal carbon dioxide tension. Thus, the remaining unknown parameters are ρ and PminO_{2}.

By definition PminO_{2} must lie between 0 mmHg and the oxygen tension of the capillary bed, PcapO_{2}. In vivo studies suggest the oxygen tension at the brain mitochondria is low in healthy brain ^{17, 30}, with an average value of approximately 8.5 mmHg. In this work we use a gamma distribution to model PmitO_{2}, maintaining the same median value but allowing variation between zero and PcapO_{2}. Thus, including the entire range of feasible values for the mitochondrial oxygen tension in the healthy and diseased brain.

The in-vivo variation in ρ has not been studied directly, however, we can gain some insight by considering the variation in vascular volume fractions with ageing and disease. For example, studies of capillary density in the human brain report age related decreases of approximately 16%, while studies in rat suggest that such vascular changes are mirrored by alterations in venular blood volume, with an approximate 20% decrease in venular blood observed across the life span of rats ^{31}. Large reductions in vascular density have been observed in aging patients with neurodegeneration. For example, Buee et al. ^{32} reported a 38% decrease in vessel density in elderly AD patients compared to a healthy 49-year-old. However, the fractional change in the capillary density is reported to be of a similar magnitude, suggesting minimal alterations in ρ ^{33, 34}.

MRI studies of stroke using vessel size imaging (VSI) report an increase in the mean vessel size and a reduction in vessel density. Although histological studies confirm a significant reduction in vessel density, they do not report any change in the mean vessel size ^{35}. The in-vivo observation of an increase in vessel size is likely to be due to vasodilation of arterioles, and thus the capillary fraction of the BOLD sensitive vasculature is likely to remain unchanged even in stroke.

Taken together these findings suggest there is little variation in the relationship between capillary blood volume and the BOLD sensitive blood volume in a wide array of physiological and pathological conditions including ageing, neurodegeneration, and stroke. In our modelling we chose the range of ρ to be 2 to 3.33, giving a capillary blood volume of 20 to 40% of total blood volume, when the arterial contribution is assumed to be 20 to 30% ^{36}. This range was chosen to offer reasonable physiological variation around a typical capillary volume of 33% ^{37, 38}.

The scaling factor, A, in the BOLD equation is similar to β, in that they both capture information related to vessel geometry and field strength ^{39}. While a fixed value is normally given to β, the value of A is normally left undefined. However, we are able to approximate A if we assume that the majority of the BOLD signal arises from the extravascular space around veins and This simplification follows that used by Blockley et al. who used R_{2}′ as a surrogate for M in calibrated fMRI studies. ^{16, 40}. Blockley et al demonstrate that R_{2}′ captures most of the BOLD signal variation and predict a tight correlation between R_{2}′ and M. Thus, with reference to equation 1, a reasonable approximation of A can be obtained from in-vivo measurements of R_{2}′, CBV_{v}, and SvO_{2}. The cortical R_{2}′ is approximately is 3 s^{−1} at 3T ^{41}, while the BOLD sensitive blood volume has been estimated to be between 1.75% ^{6} and 3.6% ^{42}. Assuming an average [Hb] of 14 g/dL, an SvO_{2} of 0.6, and taking the mean CBV_{v} of 2.5% we obtain a value of 14 s^{−1}g-^{β}dL^{β} for A at 3T.

By taking this estimate of A as the mean value and including variance to capture the physiological uncertainty, we are able to perform Monte-Carlo simulations to explore the relationship between M, CBF_{0} and OEF_{0}. We performed these simulations with both the standard definition of M and the newly proposed definition (equation 5). Table 1 outlines the physiological parameters and values used in the Monte-Carlo simulations.

When oxygen diffusivity is not included in the definition of M, but the physiology is constrained as per table 1, we are able to fit a Lowess surface (quadratic, span 25) to OEF_{0} (from M and CBF_{0}) with an RMSE of 0.086 (R^{2} of 0.70) figure 1a. If oxygen diffusivity is included in the definition of M, figure 1b, we are able to fit a surface with an RMSE of 0.043 (R^{2} of 0.92), significantly reducing the error in OEF_{0}. Thus, we are able to predict the modelled OEF_{0} with a RMSE error of approximately 10%, using only measures of M and CBF_{0} (given a fixed [Hb]).

Because the Monte-Carlo simulations in figure 1 have a fixed [Hb] we cannot generalize the relationship across Hb concentrations. However, if a surface is fit for CMRO_{2,0} rather than OEF_{0}, then [Hb] does not have to be considered as a separate parameter (at least under normoxic conditions). Figure 2 shows how a simple polynomial surface (order 2,3) can be fitted to M and CBF_{0} to estimate CMRO_{2,0} across a wide range of Hb values (10 to 18 g/dL) with a RMSE of 22.5 μmol/100g/min (R^{2} = 0.94).

The inverse relationship, mapping CBF_{0} and CMRO_{2,0} to M, can also be estimated with a surface. Thus, allowing the maximum BOLD signal to be directly estimated from the physiological parameters. Equation 6 shows a rational equation that allows for such a mapping (RMSE = 0.04, R^{2} = 0.85).

The modelling suggests that by fitting for M and CBF_{0}, we can make estimates of CMRO_{2,0} with a low degree of uncertainty, across a wide range of physiology, and using only a single hypercapnic challenge. However, as highlighted by equation 3b, the BOLD response to breath-hold challenges has a dependence on CMRO_{2,0}. Therefore, we cannot estimate CMRO_{2,0} in this step like manner (without making further approximations) and must instead make simultaneous estimates of M and CMRO_{2,0}, for example, with numerical methods. Here we employ a gradient boosted decision tree algorithm, LightGBM, which is trained on simulated data to learn the mapping between MRI data and CMRO_{2,0}. Thus, allowing rapid estimation of CMRO_{2,0} directly from the MRI data. The details of the machine learning analysis method and training are given in section 3.

## 3 Machine Learning Methodology

The machine learning regressors employed in this work are trained on simulated data. The method and its implementation are similar to our recently presented work on robust parameter estimation with dc-fMRI ^{19}. In the present implementation the BOLD signal model (equations 3b and 5) and simplified pCASL kinetic arterial spin labelling signal model ^{44} were used to generate artificial MRI time series to match the in-vivo acquisition and breath-holding protocols.

The relaxation rate of arterial blood, R_{1}, was modelled as in shown in equation 7 ^{45, 46} and the P50 for Hb saturation was estimated from the partial pressure of resting end-tidal CO_{2}, as described in ^{45}.

An acquisition length of 8 minutes 45 seconds was simulated with a TR of 4.4 seconds (119 volumes) and a BOLD echo time of 30ms. The breath-holding paradigm is detailed in figure 3, and table 2 lists the values and distributions of all MRI and physiological parameters used in the data simulations. Breath-holds were modelled using a boxcar function convolved with a gamma density function. Both hypercapnic and hypoxic variation was modelled, each with separate rise and fall times, taking into account the simultaneous changes in PaCO_{2} and PaO_{2} during breath-holding ^{26}. Variation in local arrival times was modelled with a bulk delay of up to 3 TR’s (13.2 seconds).

Simulated MRI data was high pass filtered (300 second cut-off for the BOLD data only), and Fourier transformed to create frequency domain ASL and BOLD data. A feature vector was constructed by concatenating physiological and sequence parameters; [Hb], CaO_{2,0}, post-label delay, and the first 15 points of the frequency domain data (excluding the DC BOLD component). The implementation in this work used only the MRI magnitude data (no phase data). This is a change to the previous implementation where both magnitude and phase data were used ^{19}. The motivation for this change is to reduce the sensitivity to hemodynamic delays and is practicable due to the periodic nature of the breath-hold stimulus. We used LightGBM ^{20}, an efficient gradient boosted decision tree, to learn the mapping between the MRI data and the two target parameters, CBF_{0} and CMRO_{2,0}. LightGBM is a state-of-the-art gradient boosted decision tree algorithm that is computationally efficient and is suitable for training on big data sets. A summary of the machine-learning pipeline is presented in figure 4.

50,000 simulations were used for training an individual regressor to predict CBF_{0} and CMRO_{2,0}. Training of each network took less than 10 seconds on a 2013 MacBook Pro with 16GB of memory. We used the Python implementation of LightGBM (version 3.0.0.99) for training with the maximum number of leaves per tree set to 1,000; all other parameters were default. To assess the performance of the LightGBM implementation 10,000 additional simulations were performed using the distribution of parameters detailed in table 2, but with OEF_{0} limited to 0.15 to 0.65.

## 4 In-vivo experiments

Fifteen healthy volunteers (4 males, mean age 27.4 ± 6.2 years) were recruited to the study. The local ethics committee approved the study and written informed consent was obtained from each participant. Blood samples were drawn via a finger prick prior to scanning and were analysed with the HemoCue Hb 301 System (HemoCue, Ängelholm, Sweden) to calculate the systemic [Hb] value for each participant. All MRI data were acquired using a Siemens MAGNETOM Prisma (Siemens Healthcare GmbH, Erlangen) 3T clinical scanner with a 32-channel receiver head coil (Siemens Healthcare GmbH, Erlangen). An 8 minute 45 second dual-excitation pseudo-continuous arterial spin labelling (pCASL) and BOLD-weighted acquisition was acquired during repeated breath-holding (see figure 3 for experimental protocol). End-tidal gas monitoring was performed throughout the acquisition via a nasal cannula using a rapidly responding gas analyser (PowerLab®, ADInstruments, Sydney, Australia). All volunteers practised the breath-hold task to ensure they understood the visual instructions, prior to the MRI session.

The parameters for the in-house pCASL sequence ^{45} were as follows: post-labelling delay and label duration 1.5 seconds, EPI readout with GRAPPA acceleration (factor = 3), TE_{1} = 10ms, TE_{2} = 30ms, TR = 4.4 seconds, 3.4 x 3.4 mm in-plane resolution, and 15 (7mm) slices with 20% slice gap. A 3D magnetization-prepared rapid gradient-echo (MPRAGE) sequence was acquired for image registration purposes (1mm slice thickness, 1.14 × 1.14 in-plane image resolution, TR/TE = 2100/3.2 ms). FAST ^{47} was used for tissue segmentation to create high-resolution grey matter partial volume estimates.

Post processing of the pCASL data included motion correction of the time series with the FSL tool MCFLIRT ^{48} followed by surround subtraction and surround averaging of TE_{1} and TE_{2} time series to create perfusion-weighted and BOLD-weighted time series respectively. After temporal smoothing, the BOLD time series was high pass filtered with a 300 second cut-off using the filter implementation in FSL. Each time series was then Fourier transformed, and the first 15 points of the magnitude data were included in a feature vector used for parameter estimation. The DC component from the BOLD data was excluded from the feature vector such that only 14 BOLD data points were included. The feature vector was completed with the inclusion of CaO_{2,0}, [Hb], and the post-labelling delay (calculated per slice). CBF_{0} and CMRO_{2,0} parameter estimation was performed using the LightGBM estimator trained on simulated data.

Subjects’ resting perfusion estimates were used to register the low-resolution functional data to high-resolution grey matter partial volume estimates. An inverse transform was used to create partial volume estimates in functional space. A grey matter partial volume estimate of 0.5 was taken as the threshold to create grey matter masks for statistical reporting of grey matter averages.

## 5 Results - Simulations

Cross-validation of the regression models show that CBF_{0} estimates have an R^{2} of 0.99 and a RMSE of 0.3 ml/100g/min, while CMRO_{2,0} estimates have an R^{2} of 0.94 and a RMSE of 22.9 μmol/100g/min. This is a similar level of uncertainty observed in section 2, when using known M and CBF_{0} values, demonstrating the ability of the regressor to estimate CMRO_{2,0} directly from the simulated MRI data. In addition to quantifying the RMSE for individual voxels it is equally important to assess the bias in parameter estimates and their sensitivity to variation in cerebral physiology. Figure 5 shows the percentage error in CMRO_{2,0} estimates for key modelling parameters. The error plots typically have a bias of less than 10%, with maximum values approaching 20% at the extreme ranges of parameter values.

The modelling parameter with the largest bias is PmitO_{2}. However, if we plot a surface of the mean percentage CMRO_{2,0} error against MTT and PmitO_{2}, figure 6, we see that the most significant bias only occurs when MTT is long and PmitO_{2} is high. Although this is potentially problematic in certain pathologies, it is unlikely to present a practical limitation, as ASL imaging at 3T is unlikely to provide useful data when transit times are very long.

Taken together, these results demonstrate that the implemented ML algorithm is able to accurately capture the modelled relationship between the MRI data and the physiological parameters across a wide distribution of physiological parameters. The lack of significant biases suggests the method should be sensitive to metabolic variation in both healthy and pathological tissue when applied in-vivo.

## 5 Results - In Vivo

All subjects were able to follow the experimental instructions and successfully performed the respiratory task. We were unable to obtain a blood sample from one subject out of fifteen scanned. This subject was excluded from further analysis, as they did not have a measured [Hb]. The average grey matter temporal signal (tSNR) to noise was 2.8 for the ASL acquisitions and 90 for the BOLD acquisitions.

Table 3 provides a summary of average grey matter results each for subject. M values were calculated voxelwise using equation 6 to find an M value consistent with the estimated CBF_{0} and CMRO_{2,0}. Matching parameter maps from each subject are displayed in figure 7.

Figure 7 shows example CBF_{0}, OEF_{0}, CMRO_{2,0}, and M maps for each subject. As can be seen from the figure, and consistent with expected physiology, regions of high perfusion are matched with areas of high metabolic oxygen consumption, while the oxygen extraction fraction has little variation throughout the grey matter. However, OEF_{0} estimates in pure white matter voxels are artificially high. This is likely to be a consequence of the long arrival time of blood in white matter leading to a reduced perfusion signal. In accordance with expectations, the M maps contain higher values in grey matter compared to white matter (due to increased blood volume). However, due to the overestimation of OEF_{0} in white matter M is also overestimated here, reducing the expected grey matter to white matter contrast.

A regression of grey matter OEF percentage against [Hb] and CBF produced the following model (p < 0.05).

The parameter fits are significant for both CBF and [Hb] (p<0.05) and are in close agreement with the results reported by Ibaraki et al 2010 ^{49} and with our previous work using dual-calibrated fMRI. This result suggests that the magnitude of parameter estimates made with the proposed method are appropriate, and that the method is able to detect underlying variations in physiology and metabolism, as predicted by the modelling.

## 5 Discussion

This work presents a practical and pragmatic method for mapping the rate of cerebral grey matter metabolic oxygen consumption with MRI. The approach is straightforward to apply and should be widely applicable to research as well as clinical studies. Although the approach requires validation with proven measures of cerebral oxygen metabolism, the agreement between in-vivo results and previously observed physiological relationships is encouraging and supports the validity of the method.

We have chosen to demonstrate the method using a breath-holding protocol with simultaneously acquired BOLD and ASL data. This method provides rapid measurement of both M and CBF_{0}, which are required for the calculation of oxygen extraction and CMRO_{2}. We have also implemented a machine-learning analysis approach for rapid parameter estimation. However, the method can in principle be applied to any data that can provide a robust measurement of M and an estimate of resting CBF. To facilitate such analysis, we have provided a simple mapping from M and CBF_{0} to CMRO_{2,0}.

We show robust parameter estimates in a healthy young cohort based on a breath-hold of 20 seconds, repeated 10 times, and interleaved with paced breathing. Many breath-hold repeats are needed to ensure sufficient signal to noise when modelling ASL data with a long effective TR. Despite some extra challenges that may arise for clinical applications, this breath-hold method is non-invasive and does not require specialist respiratory equipment for gas inhalation, therefore it has potential for wider applicability in research and clinical settings. Breath holding tasks during MRI have been used successfully in many clinical populations ^{50–56}, mostly during BOLD fMRI to map cerebrovascular reactivity. Task cues and practice sessions may need to be tailored to each clinical population to ensure adequate compliance.

The limitations of the method are mostly shared with the dual-calibrated methodology, see ^{57} for a review. For example, for the method to be viable there must be a local increase in CBF with breath-holding, i.e., there must be a local vascular reserve. This condition may not be met in diseases such ischemic stroke, where arterial vessels are maximally dilated in an attempt to maintain local perfusion. Additionally, the method is expected to be vulnerable to large changes in the ratio of the capillary to deoxyhaemoglobin sensitive blood volumes, i.e. beyond the typical variation in vascular compartments referenced in this work. Although large changes in this ratio appear unlikely in many brain diseases, we might expect such vascular remodelling to occur in advanced or aggressive brain tumours ^{58}. Therefore, we would encourage caution in interpreting the results in such applications. Nevertheless, the proposed method offers a simple means of mapping cerebral oxygen metabolism with MRI and has the potential to be a useful tool for both neuroscience research and clinical imaging.

## Data Availability Statement

The python code for the data analysis pipeline proposed in this article is freely available https://zenodo.org/badge/latestdoi/350795130. We do not have ethical consent to make the in-vivo datasets acquired for this study publicly available.

## Acknowledgements

We would like to thank the UK Engineering and Physical Sciences Research Council (EP/S025901/1) and Wellcome for supporting this work: Wellcome Strategic Award, Multi-scale and multi-modal assessment of coupling in the healthy and diseased brain, grant reference 104943/Z/14/Z. MG thanks the Wellcome Trust for its generous support via a Sir Henry Dale Fellowship (220575/Z/20/Z).