Dual-calibrated fMRI measurement of absolute cerebral metabolic rate of oxygen consumption and effective oxygen diffusivity.

Dual-calibrated fMRI is a multi-parametric technique that allows for the quantification of the resting oxygen extraction fraction (OEF), the absolute rate of cerebral metabolic oxygen consumption (CMRO2), cerebral vascular reactivity (CVR) and baseline perfusion (CBF). It combines measurements of arterial spin labelling (ASL) and blood oxygenation level dependent (BOLD) signal changes during hypercapnic and hyperoxic gas challenges. Here we propose an extension to this methodology that permits the simultaneous quantification of the effective oxygen diffusivity of the capillary network (DC). The effective oxygen diffusivity has the scope to be an informative biomarker and useful adjunct to CMRO2, potentially providing a non-invasive metric of microvascular health, which is known to be disturbed in a range of neurological diseases. We demonstrate the new method in a cohort of healthy volunteers (n = 19) both at rest and during visual stimulation. The effective oxygen diffusivity was found to be highly correlated with CMRO2 during rest and activation, consistent with previous PET observations of a strong correlation between metabolic oxygen demand and effective diffusivity. The increase in effective diffusivity during functional activation was found to be consistent with previously reported increases in capillary blood volume, supporting the notion that measured oxygen diffusivity is sensitive to microvascular physiology.

Introduction

Calibrated fMRI measurement of absolute cerebral rate of oxygen metabolism (CMRO2) offers a non-invasive method of mapping oxygen consumption in the brain (Bulte et al., 2012; Gauthier et al., 2012; Wise et al., 2013), providing quantitative estimates of a critical physiological function. However, the method does not directly consider the transport of oxygen into the tissue, which is principally constrained by cerebral blood flow (CBF) and the effective oxygen diffusivity of the capillaries (Buxton and Frank, 1997; Gjedde et al., 1999; Hayashi et al., 2003; Hyder et al., 1998; Mintun et al., 2001; Vafaee and Gjedde, 2000; Valabregue et al., 2003; Zheng et al., 2002). Effective oxygen diffusivity summarises the practical ability of the capillary network for oxygen diffusion into the tissue and limits the speed of oxygen transport out of the microvasculature. One of the primary determinants of the effective oxygen diffusivity is the capillary density (Gjedde et al., 1999), which is known to be associated with mitochondrial density (Hoppeler and Kayar, 1988) and metabolic demand (Harrison et al., 2002). Thus, brain regions with a high resting CMRO2 are found to be co-localised with regions of high capillary density (Sette et al., 1989). However, the effective diffusivity does not appear to be a fixed property of the tissue and may play a crucial role in neurovascular coupling, with oxygen diffusivity being observed to parallel increases in demand and compensate for reductions in oxygen delivery (Hayashi et al., 2003, 2004; Hyder et al., 1998; Vafaee and Gjedde, 2000).

Compartmental models of oxygen exchange between the capillaries and tissue offer a means of estimating the effective oxygen diffusivity from observations of blood flow and oxygen extraction. The model proposed by (Hyder et al., 1998) suggests a need for the effective diffusivity to increase during functional hyperaemia in order to meet the metabolic demands of neural activation. Based on a meta-analysis of CBF and CMRO2 measurements from a variety of modalities Hyder et al. proposed a linear coupling between flow and effective diffusivity to account for this apparently coupled behaviour. Evidence for this linear relationship between the effective oxygen diffusivity and CBF was demonstrated with a combined MRI and MRS approach in rat (Hyder et al., 2000). However, PET experiments conducted by (Vafaee and Gjedde, 2000) demonstrate a need for the oxygen diffusivity to adapt to the current metabolic demand, with alterations in the effective diffusivity appearing to be made independently from cerebral blood flow. Alternatively, more recent analysis presented by (Buxton, 2010) demonstrates that metabolic oxygen demand could be met if there was fixed but significant oxygen tissue content, without the need for adjustment to the oxygen diffusivity. The exact mechanism responsible for any such adaptation to metabolic demand is unclear. However, a plausible candidate for the modulation of the effective diffusivity is via pericyte control of capillary dilation, either through a direct increase in the capillary blood volume, or via a homogenisation of flow heterogeneity (Jespersen and Ostergaard, 2012). Thus, measurement of the effective diffusivity may in fact provide a non-invasive probe to investigate the health and action of capillary pericytes, whose function is known to be degraded in multiple neurodegenerative diseases and stroke (Winkler et al., 2013; Yemisci et al., 2009; Zlokovic, 2011).

In the work presented here we use a compartmental model of oxygen exchange to model the relationship between blood flow, effective diffusivity and oxygen extraction. The model is included within a dual-calibrated fMRI estimation framework (Germuska et al., 2016) to enable simultaneous estimates of the resting blood flow, oxygen extraction fraction (OEF), effective diffusivity, and CMRO2. The aim of this study was to examine the coupling between consumption (CMRO2) and diffusivity at rest, and in response to neural activation (a visual checkerboard task) using the newly proposed method. Our first hypothesis was that, due to the tight functional-structural coupling between capillary density and resting metabolism, there would be a strong correlation between the basal CMRO2 and effective diffusivity. Secondly, we hypothesised that the increased metabolic demand due to the visual task would result in a parallel increase in effective diffusivity, whose magnitude should be consistent with published recordings of functional capillary recruitment (Hall et al., 2014).

Methods

Compartmental modelling

The compartmental model of oxygen exchange is based on the model of (Hayashi et al., 2003). As shown in Fig. 1, the model contains a single capillary compartment with area A and length L, which exchanges oxygen with a cylindrical volume of tissue. The capillary has two compartments, a haemoglobin compartment (with oxygen content CB) and a plasma compartment (with oxygen partial pressure P). The oxygen in plasma is assumed to be in equilibrium with the oxygen bound to haemoglobin as described by the Hill equation belowwhere P is the partial pressure of oxygen in the plasma, P50 is the oxygen partial pressure at half saturation, [Hb] is the haemoglobin concentration (g/ml), φ is the oxygen binding capacity for Hb (1.34 ml/g), h is the Hill coefficient (2.8), and CB (oxygen bound to haemoglobin) is equal to total capillary oxygen content, Ct, if the contribution of plasma oxygen is neglected.

Fig. 1

Schematic of the simple compartmental model for oxygen exchange between capillary blood and brain tissue.

As blood travels along the capillary, the oxygen exchanges between an infinitesimally thin element of blood plasma and a well-stirred oxygen compartment some fixed distance from the capillary (with partial pressure Pm). The permeability of the capillary endothelium and brain tissue are combined into a single effective permeability, k. This interpretation of the model is a departure from the model presented by (Hayashi et al., 2003), who assumed a uniform partial pressure of oxygen in the radial slice of plasma and zero partial pressure of oxygen on the tissue side of the capillary-tissue interface, thus, localising the oxygen transport to within the capillary endothelium. However, in-vivo measurements suggest that the capillary wall does not present a significant barrier to oxygen diffusion (Duling et al., 1979), which is instead provided by the tissue (Hudetz, 1999). Thus, as per (Hyder et al., 1998), we combine both the capillary wall and the surrounding brain tissue into a single interface between the plasma and a well-stirred pool at the end of the diffusion path, which is presumably within or surrounding the mitochondria.

From our compartmental model we can define the differential equation describing the loss of oxygen from within a capillary aswhere t is time, x is the fractional distance along the capillary (0,1), k is the effective permeability (mL/mmHg/mL/min), P is the oxygen partial pressure in the plasma (mmHg), and Pm is the oxygen partial pressure at the mitochondria (mmHg).

Following (Zheng et al., 2002) and (Hayashi et al., 2003)and when blood flow is constantwhere CBF is cerebral blood flow in ml/100 g/min, A is the cross-sectional area of the capillary, L is the length, and V is the volume in ml/100g. Thus, by substituting equations (3), (4) into equation (2) we obtain

If we assume that there is minimal partial pressure of oxygen at the mitochondria (Gjedde et al., 1999, 2005; Herman et al., 2006), i.e. Pm ≈ 0, then the product k×V is the effective oxygen diffusivity, DC (ml/100g/mmHg/min). By substituting equation (1) into equation (5), we can express the differential equation of oxygen loss in terms of capillary oxygen content, resting CBF, and the effective oxygen diffusivity aswhich is the equivalent to equation (6) in (Hayashi et al., 2003), except that we are assuming negligible oxygen tension at the mitochondria rather than zero average oxygen tension in the tissue.

At the macroscopic level we consider a volume of tissue to contain a collection of identical capillaries arranged such that Pm can be considered uniform (note this does not pre-suppose any particular structural configuration). For simplicity we also assume that all other parameters are identical across the capillaries, such that there is no flow heterogeneity or variation in haemoglobin concentration [Hb]. Thus, the modelled oxygen diffusivity represents a combination of vascular parameters including capillary blood volume, flow heterogeneity, and any underlying variation in Pm.

Equation (6) was solved numerically (using MATLAB's ordinary differential equation solver (Mathworks, MA, USA.)) for different combinations of Dc, CBF, P50, and [Hb] to create a lookup table of results that could be used to fit in-vivo data. The oxygen extraction fraction was calculated by evaluating (Ct(x)|x=0 - Ct(x)|x=1)/Ct(x)|x=0 for each combination of parameters, where Ct(x)|x=0 is the oxygen content at the arterial end of the capillary (CaO2), assumed to be 0.95 of the maximum (Jespersen and Ostergaard, 2012), and Ct(x)|x=1 is the oxygen content at the venous end of the capillary (CvO2).

Calibrated fMRI signal modelling

Quantification of the oxygen extraction fraction and resting blood flow (from which CMRO2 is calculated) is performed using the dual-calibrated fMRI method (Bulte et al., 2012; Gauthier et al., 2012; Wise et al., 2013) within a forward modelling framework (Germuska et al., 2016). The method is based upon the isometabolic alteration of flow and venous oxygenation using hypercapnic and hyperoxic respiratory modulations. Here we utilise the simplified calibration model (Merola et al., 2016), where the change in BOLD signal is defined by equation (7).where S is the MR signal magnitude, TE is the echo time of the acquisition, κ is a composite calibration parameter that represents the combination of the venous-weighted blood volume and water diffusion effects, [dHb] is the deoxyhaemoglobin concentration and is equal to [Hb]×(1-SvO2), θ (assigned a value of 0.06) is an empirical parameter combining contributions from venous blood volume changes during hypercapnic hyperaemia and extra-vascular water diffusion effects around the microvasculature, and the subscript 0 represents the a parameter's baseline value.

The deoxyhaemoglobin ratio is modelled as shown in equation (8) (Wise et al., 2013), and as before, OEF = (CaO2 – CvO2)/CaO2.

The arterial spin labelling (ASL) sequence used for the calibrated acquisition uses a pCASL labelling scheme with pre-saturation and background suppression (Okell et al., 2013), and a dual-excitation (DEXI) EPI readout (Schmithorst et al., 2014). As such it differs from the dual-echo PASL acquisition previously employed in the forward modelling framework (Germuska et al., 2016), and the methods have been adapted to reflect this. In the current implementation, BOLD contamination is removed from TE1 via surround subtraction, and ASL contamination is removed from TE2 via surround averaging. Thus, only the BOLD model (equations (7), (8) is used to estimate TE2 data, while TE1 time courses are estimated according to the simplified pCASL kinetic model (Alsop et al., 2015), equation (9).where ΔS is the tag/control difference, α is the tagging inversion efficiency (0.85), αinv is a scaling factor to account for the reduction in tagging efficiency due to background suppression (0.88)(Mutsaerts et al., 2014; Shin et al., 2011), T1,blood is the longitudinal relaxation time of arterial blood, M0 is the equilibrium magnetisation, λ is the brain/blood partition coefficient (0.9), τ is the tagging duration, and PLD is the post labelling delay. See Table 1 for a summary of the parameters used in the modelling of ASL, BOLD and oxygen exchange.

Table 1

Abbreviations for variables and techniques used in the modelling and analysis.

Variable/abbreviation	Expression (units)
OEF	Oxygen Extraction Faction (dimensionless 0–1)
CMRO2	Cerebral Metabolic Rate of Oxygen consumption (μmol/100 g/min)
CBF	Cerebral Blood Flow (ml/100g/mmHg/min)
CVR	Cerebral Vascular Reactivity (% CFB change/mmHg Co2)
Dc	Effective oxygen diffusivity of the capillary network (ml/100g/mmHg/min)
P	Oxygen tension in capillary plasma(mmHg)
P50	Oxygen tension at which haemoglobin is 50% saturated (mmHg)
Pm	Oxygen tension at the mitochondria (mmHg)
[Hb]	Haemoglobin concentration (g/ml)
CB	Oxygen content bound to haemoglobin (ml/ml)
Ct	Total capillary oxygen content (ml/ml)
CaO2	Oxygen content at the arterial end of the capillary network (ml/ml)
CvO2	Oxygen content at the venous end of the capillary network (ml/ml)
SaO2	Arterial oxygen saturation (dimensionless 0–1)
SvO2	Venous oxygen saturation (dimensionless 0–1)
PaO2	Arterial oxygen tension (mmHg)
PaCO2	Arterial carbon dioxide tension (mmHg)
PETO2	End-tidal oxygen tension (mmHg)
PET-CO2	End-tidal carbon dioxide tension (mmHg)
ϕ	Oxygen binding capacity of haemoglobin (1.34 ml/g)
h	Hill coefficient (2.8)
k	Effective permeability of capillary endothelium and brain tissue (ml/mmHg/ml/min)
ε	Oxygen plasma solubility (0.0031 ml/mmHg/dl)
BOLD	Blood oxygenation level dependent MRI signal
ASL	Arterial spin labelling
TE	Echo time of MRI acquisition (ms)
κ	BOLD calibration parameter including venous-weighted blood volume and water diffusion effects
[dHb]	Deoxyhaemoglobin concentration (g/ml)
θ	Effective hypercapnic venous flow-volume coupling constant (0.06)
T1,blood	Longitudinal relaxation time of arterial blood (s)
R1,blood	Longitudinal relaxation rate of arterial blood (s−1)
M0	MRI signal equilibrium magnetisation (dimensionless)
λ	Brain/blood partition coefficient (dimensionless, 0.9)
τ	Arterial spin labelling tagging duration (s)
PLD	Arterial spin labelling post labelling delay (s)

Data acquisition

Nineteen healthy volunteers (13 males, mean age 31.9 ± 6.5 years) were recruited to the study. Volunteers' tolerance of hypercapnic periods and breathing through a face-mask was tested with a non-MRI session prior to MRI scanning. The study was approved by the local ethics committee. Written informed consent was obtained from each participant. All 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). During each scanning session an 18-min dual-calibrated fMRI scan was acquired with interleaved periods of hypercapnia, hyperoxia and medical air being delivered to the subjects according to the protocol previously proposed by our lab (Germuska et al., 2016). End-tidal gases, PETCO2 and PETO2, were sampled from the volunteer's facemask using a rapidly responding gas analyzer (AEI Technologies, Pittsburgh, PA, USA), see Fig. 2 for a summary of end-tidal recordings and timings of the gas paradigm.

Fig. 2

Mean (solid line) and standard deviation (shaded area) of end-tidal recordings from all subjects included in analysis (n = 16). Absolute value of end-tidal oxygen partial pressure (red) and relative change in end-tidal carbon dioxide partial pressure (mmHg).

All calibrated fMRI data were acquired using a prototype pCASL acquisition using pre-saturation and background suppression (Okell et al., 2013) and a dual-excitation (DEXI) readout (Schmithorst et al., 2014), see Fig. 3 for a sequence timing diagram. The labelling duration and PLD were both set to 1.5s, GRAPPA acceleration (factor = 3) was used with TE1 = 10 ms and TE2 = 30 ms. An effective TR (total repetition time including labelling scheme and both readout periods) of 4.4 s was used to acquire 15 slices, in-plane resolution 3.4 × 3.4 mm and slice thickness 7 mm with a 20% slice gap. A calibration (M0) image was acquired for ASL quantification with pCASL and background suppression switched off, with TR of 6 s, and TE = 10 ms.

Fig. 3

Pulse sequence timing diagram for dual-excitation pseudo-continuous arterial spin labelling (DEXI-pCASL) acquisition. Sequence timings are in ms to the nearest 5 ms.

For a subset of volunteers (n = 7, 3 males, mean age 37.4 ± 6.7 years) an additional 8-min black and white visual checkerboard task (reversing at a frequency of 2Hz, alternating between 30 s rest and 30 s stimulus) was performed during pCASL DEXI data acquisition. In each volunteer a T1-weighted structural image was acquired for grey matter segmentation and to aid registration to standard (MNI) space.

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. The partial pressures of end-tidal gas concentrations were assumed to be in equilibrium with arterial blood, such that PaO2 = PETO2 and PaCO2 = PETCO2. Baseline PaCO2 recordings were used to estimate resting blood pH based on the Henderson-Hasselbalch equation (equation (10)), assuming HCO3- = 24 mmol/L (Gai et al., 2003).

From which the resting P50 was calculated according to the linear correlation, P50 = 221.87–26.37×pH, reported by (Gai et al., 2003). The Severinghaus equation (Severinghaus, 1979) was used to convert PaO2 recordings into SaO2 time series, which were then converted to CaO2 via equation (11).where ε is the O2 plasma solubility (0.0031 ml/mmHg/dL). The T1 of arterial blood was calculated from a linear fit to SaO2, PaO2 and R1,blood in-vivo data presented in (Pilkinton et al., 2012), equation (12).where R1,blood is the longitudinal relaxation rate of arterial blood in seconds.

Data analysis

Data were pre-processed using a combination of MATLAB code and FSL (Jenkinson et al., 2012). Motion correction was performed with the FSL MCFLIRT function and spatial smoothing (FWHM = 4.5 mm) of the BOLD data (surround average of TE2) was performed with SUSAN (Smith and Brady, 1997). ASL data (surround subtraction of TE1) and M0 acquisition were spatially smoothed using a 3D Gaussian kernel (FWHM 4.5 mm). DEXI data was registered to the structural T1 data using FSL's epi-reg tool. Following grey matter segmentation of the structural T1 image, using FAST (Zhang et al., 2001), grey matter estimates were transformed to native space and used for grey matter masking (threshold of 0.5). DEXI data was masked prior to analysis using a binarised M0 image to reduce processing time.

End-tidal traces were aligned with the DEXI data via a cross-correlation between PaCO2 and the mean grey matter ASL signal. Measured [Hb] and calculated P50 values were used to resample the initial 4D effective diffusivity lookup table to a high-resolution 2D lookup table relating CBF0 and Dc to OEF, enabling simple linear interpolation to be used during fitting. MATLAB's non-linear least squares minimisation routine (lsqnonlin) was used to simultaneously optimise voxelwise estimates of Dc, OEF0, CBF0, κ, and the cerebral vascular reactivity (CVR) by minimising the least squares difference between the acquired data and modelled ASL and BOLD timeseries (MATLAB code for pre-processing of end-tidal traces and parameter estimation is available from 10.5281/zenodo.1285862 and 10.5281/zenodo.1285845). See Fig. 4 for a flow diagram representing the forward model used in the analysis framework. In line with the forward modelling approach previously published (Germuska et al., 2016) regularisation was applied to reduce instability in fitting a non-linear model to the data (see appendix for full details). Briefly, regularisation was applied to DC and OEF in an adaptive manner to reduce the sensitivity to noise variation across voxels and subjects. The regularisation parameter for oxygen extraction fraction was assumed to be uniform, with a nominal OEF0 of 0.4. We make the assumption that oxygen diffusivity varies with capillary density (and therefore grey matter partial volume) to impose spatial variation on the diffusivity regularisation. Grey matter partial volume estimates were calculated by normalising an initial perfusion estimate by its maximum value (median value in 100 voxels with greatest signal intensity) and then multiplying by 0.15 ml/100g/mmHg/min. We use an initial perfusion estimate (rather than a segmented structural image) to estimate grey matter partial volume to avoid bias due to segmentation and registration errors. As per our previous work (Germuska et al., 2016), we used digital phantom experiments to determine the optimal level of regularisation for each parameter, OEF and DC. Additionally, we explored the influence of the SNR on the mean squared error of regularised fits to the simulated data; see the appendix for further details of the simulations.

Fig. 4

Flow diagram showing how measured physiological data and estimated parameters are combined to estimate ASL and BOLD signal time courses during parameter estimation. The forward model incorporates oxygen diffusivity modelling into a dual-calibrated fMRI framework.

Visual data was subject to the same pre-processing steps as the baseline data, additionally it was registered (FLIRT) to baseline data to account for any gross subject motion between the datasets. Percentage change in CBF and the BOLD signals were calculated using FEAT to fit the pre-processed data with the visual paradigm. CMRO2 was calculated both on a voxel-wise basis, and from a grey matter ROI, which was thresholded to include only voxels with significant BOLD and CBF activation (z-stats > 2.3). Because resting data was quantified using a simplified BOLD model (Merola et al., 2016), the standard calibration model (Davis et al., 1998) was modified as per equation (13) to calculate the visual CMRO2.

Estimates of CMRO2 and CBF during visual activation were used to calculate OEF via the Fick principle, while estimates of DC were made by inverting the look-up table and assuming constant [Hb] and P50.

Results

The mean and standard deviation of baseline PETO2 and PETCO2 were 116 ± 5 mmHg and 41.6 ± 3.6 mmHg, the hyperoxic respiratory modulation resulted in an average PETO2 of 325.2 ± 12.8 mmHg, while hypercapnia produced an average PETCO2 of 51.7 ± 3.5 mmHg. Three subjects did not return to a stable PETCO2 baseline during the medical air periods of the DEXI acquisition; this was judged to be a deviation of greater than 4 mmHg below the starting value. These subjects were excluded from further analysis. The mean increase in grey matter CBF during hypercapnia was 24.0 ± 3.7%. Group average values of the resting grey matter physiological parameters are reported in Table 2.

Table 2

Mean (±standard deviation) systemic and grey matter estimates at baseline, (n = 16).

[Hb] g/dl	P50 mmHg	CBF ml/100 g/min	OEF	CMRO2 μmol/100 g/min	DC ml/100g/mmHg/min	CVR %/mmHg
14.3 ± 1.5	27.1 ± 0.1	55.6 ± 6.3	0.38 ± 0.04	157.4 ± 12.3	0.092 ± 0.009	2.4 ± 0.4

The mean grey matter value of the effective diffusivity was 0.092 ± 0.01 ml/100g/mmHg/min, or 3.62 ± 0.39 μmol/100g/mmHg/min, which is in good agreement with the PET measurements made by (Ibaraki et al., 2010) and (Vafaee and Gjedde, 2000) who report value of 3.38 and 4.09 respectively. Baseline values of CMRO2 are strongly correlated with the effective oxygen diffusivity (R2 = 0.81, p < 0.01), across participants, as shown in Fig. 5a. This coupling between CMRO2 and oxygen diffusivity is also demonstrated within subjects, with a clear spatial similarity between the baseline parameter maps. Example resting parameter maps (CBF0, CMRO2,0 and Dc,0) for an individual subject are shown in Fig. 6. The similarity between the parameter maps is clear, with the maps of effective diffusivity closely following CMRO2.

Fig. 5

Scatter plots of whole brain grey matter parameter estimates at baseline for each subject (n = 16). Top panel (A) demonstrates a strong correlation between baseline metabolic oxygen consumption and effective oxygen diffusivity. Bottom panel (B) shows a strong negative correlation between oxygen extraction and oxygen delivery, such that delivery is elevated when OEF is low.

Fig. 6

Example baseline (CBF0, CMRO2,0 and Dc,0) parameter maps for an individual subject. The spatial similarity between oxygen diffusivity and the basal rate of oxygen metabolism is evidence of a strong structural-functional coupling between the two parameters in the basal state.

While the effective diffusivity appears to be strongly coupled to resting demand, we observed a strong negative correlation between baseline oxygen delivery (CaO2×CBF) and oxygen extraction (R2 = 0.57, p < 0.01) (Fig. 5b). This is consistent with previous findings of a strong negative correlation between OEF0 and CBF0 (Germuska et al., 2016; Lu et al., 2008) when [Hb] was assumed to be constant. By incorporating a measured [Hb] in the calculation of CaO2, the correlation with oxygen delivery is revealed; demonstrating the action of the cerebrovascular system to maintain a tightly controlled resting metabolic rate of oxygen consumption across participants.

To investigate the observed linear relationship between CMRO2,0 and DC,0, we solved the inverse model with fixed [Hb] (14.3 g/dl), P50 (27.1 mmHg), and CaO2,0 (0.189 ml O2/ml blood). By doing so we are able to examine the modelled co-variance between CMRO2,0 and CBF0 in isolation of other physiological factors, see Fig. 7. The analysis suggests a small positive correlation between CBF0 and CMRO2,0 is required to produce a strictly linear relationship between CMRO2,0 and DC,0, which agrees well with the in-vivo data (p < 0.05 for a paired t-test comparison between the predicted CBF0 and the normalised CBF0 estimates (CaO2,0 x (CBF0/0.189)). However, if the effective oxygen diffusivity is considered to be a constant physiological parameter, as is sometimes assumed (Vafaee et al., 2012), then a significant non-linear (exponential) relationship is predicted between CBF0 and CMRO2,0. Finally, we explored a hypothetical scenario in which the inverse of the correlation observed in this study exists between CMRO2,0 and DC,0. This scenario requires a much larger increase in CBF0 with CMRO2,0, but it is still mathematically and physiologically feasible to produce such a relationship. This analysis highlights the fact that even though DC,0 and CMRO2,0 appear to be coupled (in the examined cohort of young healthy volunteers) they provide complementary information, with CMRO2 reporting on the rate of oxygen consumption and DC reporting on the ability of the capillary network to supply oxygen to the tissue/mitochondria. The coupling between these parameters is likely the result of the expected correlation between microvascular structure (diffusivity) and function (metabolism).

Fig. 7

Modelled relationship between CMRO2,0 and CBF0 for linear increase in DC,0 with CMRO2,0 (blue), constant DC,0 (orange), and linear decrease in DC,0 with CMRO2,0 (yellow). In-vivo CMRO2,0 and normalised CBF0 (CaO2,0⋅CBF0/0.189), mean grey matter values overlaid (circles).

During visual simulation cerebral blood flow increased within the defined ROI (in the primary visual cortex, as expected) by 21.4 ± 4.6%, while CMRO2 increased by 15.1 ± 3.5%, resulting in a CBF to CMRO2 coupling ratio of 1.44 ± 0.24. The effective oxygen diffusivity was found to increase by 12.5 ± 3.6%. The CMRO2 and diffusivity increases are of a similar magnitude to those observed with PET measurements using a 4Hz yellow-blue contrast-reversing checkerboard, where CMRO2 and oxygen diffusivity were found to increase by 14.9% and 9.0% respectively. Fig. 8 shows the mean (n = 7) absolute change in physiological parameters evoked by the visual stimulus overlaid onto the selected slice in standard (MNI) space. Changes in diffusivity are co-localised with changes in CMRO2, whereas changes in CBF are more widespread.

Fig. 8

Overlay of the mean absolute change in CBF, CMRO2, and effective oxygen diffusivity evoked by the visual checkerboard stimulus for n = 7 subjects.

Fig. 9 shows the correlation between the average CMRO2 and effective diffusivity within the visual ROI at baseline (crosses) and during activation (diamonds) for each participant. It is clear from the graph that the coupling between demand and diffusivity, that was previously observed in the grey matter baseline data, is preserved during activation, R2 (baseline) = 0.95, R2 (activation) = 0.96, p < 0.01 for both datasets.

Fig. 9

Summary plot of visual ROI data for CMRO2 and effective oxygen diffusivity (including resting and activation data for n = 7 subjects). A tight correlation between CMRO2 and effective diffusivity is observed both at baseline (crosses) and during activation (diamonds), indicative of a tight coupling between the effective diffusivity and oxygen demand. The dashed lines are lines of best fit (linear regression) the dotted lines connect baseline and activation data for each subject.

Under the assumption that Pm is minimal and the effective permeability is held constant during activation, the apparent change in capillary blood volume is linearly related to the effective diffusivity (see equations (5), (6)). Thus, we can calculate the apparent flow-volume coupling relationship, CBV/CBV0 = (CBF/CBF0)η, implied from the diffusivity data. Using this line of reasoning we find that a coupling constant of 0.62 ± 0.13 is required, relating the 21% flow change to an apparent 12.5% volume increase. This result agrees surprisingly well with in-vivo observations of functional capillary vasodilation. Where the 6.7% increase in capillary diameter observed by (Hall et al., 2014) was calculated to produce a 19% increase in blood flow according to Poiseuille's law. Assuming constant capillary length this observation would predict a flow-volume coupling exponent of approximately 0.75 for their data.

Discussion and conclusions

In this manuscript we have presented a novel framework for the analysis of dual-calibrated data to produce simultaneous estimates of CMRO2 and effective oxygen diffusivity. The combined analysis of these two physiological parameters has the potential to provide useful insight into the underlying metabolic and vascular responses to different brain states and disease. The method was applied at rest and in combination with a visual task. The resting data showed a tight coupling between grey matter diffusivity and the basal rate of oxygen metabolism. This result is expected in the healthy brain, as there is significant evidence of a structural link between the density of capillaries (a significant determinant of the effective diffusivity) and metabolism (Gjedde et al., 1990; Harrison et al., 2002; Sette et al., 1989).

The effective diffusivity was also found to increase during functional activation, with a 12.5% increase in diffusivity being associated with a 15.1% increase in CMRO2 and 21.4% increase in CBF. The coupling ratio between CBF and CMRO2, 1.44, is at the lower end of in-vivo observations, which typically range from 1.3 to 5 for MRI methodologies (Leithner and Royl, 2014). Thus, the change in effective diffusivity is likely to be at the higher end of the expected range (in order to meet the oxygen demands of the elevated CMRO2 in the absence of a greater increase in CBF). Possible mechanisms to provide such an increase in effective diffusivity include a direct increase in capillary blood volume (Hyder et al., 1998), a homogenisation of capillary flow heterogeneity (Jespersen and Ostergaard, 2012), a reduction in the mitochondrial oxygen tension (Gjedde et al., 2005), or a high resting tissue oxygen tension (Buxton, 2010). While each of these factors could play a role in modulating the diffusivity of the capillary network, in-vivo measurements suggest that tissue oxygenation initially increases during functional activation and then normalises to a level slightly above its resting value (Ances et al., 2001). Thus, it is unlikely that there is a significant reduction in mitochondrial oxygen tension, which would be expected to lower tissue oxygenation rather than increase it. Alternatively, a high resting mitochondrial oxygen tension would result in a small vessel-to-tissue PO2 gradient, which, as highlighted by (Hyder et al., 2000), increases the effectiveness by which CBF can increase O2 delivery to the tissue. Our model predicts that the mitochondrial oxygen tension would need to be unrealistically high, approximately 30 mmHg to explain the CBF/CMRO2 increases observed in this study. While there is uncertainty in the value of Pm in the human brain, animal studies suggest that it is between 0.1 and 10 mmHg (Herman et al., 2006), thus a significant resting mitochondrial oxygen tension is unlikely. An alternative explanation explored in this paper is that there is a direct increase in capillary blood volume, potentially mediated by capillary pericytes; which have been demonstrated to alter capillary volume independently of arteriolar dilation (Mishra et al., 2016) and appear to play a significant role in neurovascular coupling (Kisler et al., 2017). For a given rate of perfusion, an increase in capillary volume would increase the mean transit time for blood to transverse the capillary network, producing a proportional increase in the effective diffusivity, and thus enabling greater extraction of the oxygen from the capillary bed. Our data suggest that a flow-volume coupling exponent of approximately 0.62 is required in the capillaries to provide the observed increase in effective diffusivity during visual stimulation. The implied 12.5% increase in capillary blood volume agrees well with the data and analysis presented by (Hall et al., 2014) in the mouse brain, suggesting this is a plausible explanation for the observed increase in diffusivity. Indeed more recent studies in the mouse brain suggest that such an increase in capillary blood volume is well within the range of normal physiological responses, where capillary volume has been found to increase by 10–26% during functional activation depending on the baseline diameter (Ito et al., 2017).

Although modulation of capillary blood volume appears to be sufficient to provide local control of capillary diffusivity, we cannot rule out a contribution from flow heterogeneity (Jespersen and Ostergaard, 2012), which is likely to be reduced during functional activation where smaller capillaries have been observed to dilate more than larger capillaries in rat (Kleinfeld et al., 1998; Stefanovic et al., 2008). However, the influence of flow heterogeneity is known to be dependent on the transit time distribution (Angleys et al., 2015), and confounding factors such as the heterogeneity of [Hb] have not been considered in the modelling, thus it is still unclear if this theoretical model of control is realised in-vivo.

As previously discussed, there is unlikely to be a significant resting Pm in the studied cohort, meaning the assumption of negligible Pm is unlikely to impact the results. However, this may not always be true. In theory Pm could increase in the presence of significant mitochondrial dysfunction due to the lack of oxygen uptake in the mitochondria. Currently the only direct evidence we are aware of for an increase in tissue (and therefore most likely mitochondrial) oxygen tension during mitochondrial dysfunction is from simulated dysfunction (due to a cyanide infusion) in piglets (Nielsen et al., 2013). However, there is increasing evidence of mitochondrial dysfunction in a number of neurodegenerative diseases such as Parkinson's disease (Powers et al., 2008) and Alzheimer's disease (Wang et al., 2014). Therefore, it should be highlighted that reductions in basal DC may not always correspond to a purely vascular origin and may also incorporate mitochondrial dysfunction.

In line with our previously published methods for analysis of dual-calibrated data (Germuska et al., 2016; Wise et al., 2013) we have included priors to stabilise the fitting process (see appendix). In the current implementation the priors are incorporated into parameter estimates via adaptive regularisation. Whenever priors are used to guide the fitting process there is always a trade off to be made between overfitting (not enough regularisation) and underfitting (too much regularisation). Digital phantom simulations were used to optimise the amount of regularisation and balance this trade-off. The proposed method employs regularisation on two parameters, the resting OEF and the effective diffusivity. In the case of underfitting we would expect the results to closely follow the prior, thus there would be little variation in OEF or the effective diffusivity between subjects. In contrast we find that OEF is highly correlated with resting oxygen delivery, and the effective diffusivity is tightly coupled to the resting CMRO2. Thus, it is unlikely that the results are significantly affected by underfitting. While we cannot rule out overfitting of the data, the appearance of the parameter maps is physiologically plausible and they do not suffer from significant instability within the grey matter. To further explore the influence of the proposed framework on the parameter estimates we performed simulations where DC was not estimated. In this implementation OEF was estimated directly as in (Germuska et al., 2016), and thus no regularisation was/could be placed on DC. These simulations (see appendix) showed similar levels of error to the proposed method, albeit with slightly increase RMSE in OEF estimates at low tSNR.

In conclusion, we have presented an MRI method for mapping the effective oxygen diffusivity of the capillary bed in combination with metabolic oxygen consumption. The method shows good agreement with PET literature and inferred changes in capillary blood volume are in agreement with two-photon laser microscopy measurements in animals, however, direct validation of the method is still outstanding. The proposed method is non-invasive and can be performed in a short timeframe. Previous measurements of effective oxygen diffusivity suggest it may be a valuable tool to understand the brain's response to altered oxygen supply and demand. Thus, the introduction of this method could offer a useful insight into a range of conditions and diseases with altered metabolism or vascular function.