A forward modelling approach for the estimation of oxygen extraction fraction by calibrated fMRI

.


Introduction
The rate of metabolic oxygen consumption in the brain (CMRO 2 )isa valuable biomarker in the assessment of disease and treatment effect (Heiss and Herholz, 1994;Lammertsma, 1987;Tohgi et al., 1998).Clinical studies using oxygen-15 PET have shown altered oxygen metabolism in a variety of diseases and disorders (Broich et al., 1992;Silverman and Alavi, 2005).However, such examinations cannot be routinely undertaken due to the lack of availability and the risks associated with the repeated administration of ionising radiation.Thus, there is significant interest in the development of alternative methods that can be more routinely employed to map CMRO 2 .One such emerging technique is through the calibrated measurement of the MRI blood oxygenation level dependent (BOLD) signal (Bulte et al., 2012;Gauthier et al., 2012;Germuska and Bulte, 2014;Wise et al., 2013).This approach, which we shall refer to as dual-calibrated functional MRI (dcfMRI), is based on the quantification of cerebral blood flow (CBF) and the change in venous oxygenation during periods of respiratory challenge.
The dcfMRI method requires at least two independent probes of the cerebral vasculature.Typically, hypercapnic and hyperoxic gas challenges are used, with hypercapnia producing an increase in CBF and hyperoxia directly modulating the venous oxy-haemoglobin saturation.In order to estimate the resting CMRO 2 a physiological model is used to analyse the CBF and BOLD responses to the stimuli.In previous implementations of this method the CBF and BOLD data have been analysed independently, creating intermediate variables that are then fed into a physiological model to solve for the resting oxygen extraction fraction (OEF), and thus for CMRO 2 .Approaching the problem in this stepwise manner can create an unstable solution due to the propagation of errors along the analysis pipeline.In fact solutions for OEF often do not exist with this approach due to incompatible estimates of CBF and BOLD signal changes (Gauthier et al., 2012).An alternative approach is to use numerical methods to simultaneously estimate all MR and physiological parameters in a one-step solution.
In practice this can be achieved by finding the parameter set that minimises the error between a forward model and the acquired MR data.In this approach the CBF and BOLD signals are estimated simultaneously via the physiological model, with information from the CBF signal used to infer on the BOLD signal and vice versa.Therefore this approach takes greater advantage of the available data, by combining the CBF and BOLD information.Thus, the proposed method is expected to improve upon the accuracy and repeatability of dcfMRI.In addition, all parameter estimates made with this approach can be guaranteed to exist by restricting the search space of the forward model.
In the work presented here the optimisation problem is posed as a non-linear least squares problem with L2 parameter regularization.Because the new method estimates all parameters simultaneously it becomes possible to factor out OEF from the BOLD calibration parameter (M).In addition, the forward model is made more general by replacing the fractional change in CBF with the relative cerebral vascular reactivity (CVR).These substitutions mean that the effect of regularization is not dependent on the stimulus or experimental design.Investigation of the forward model also demonstrates how the resting venous cerebral blood volume (CBV) can be estimated from the BOLD calibration parameter and the resting OEF.
Simulations using standard noise modelling techniques are performed for a cross-validation study to optimise the regularization tuning parameter (λ) and subsequently for a comparison with existing analysis methods (Bulte et al., 2012;Gauthier et al., 2012).
The tuning parameter is optimised over a wide range of baseline physiological states to minimise the bias associated with regularization.In addition to simulation studies, a comparison against existing analysis methods is made using in vivo data (acquired in 10 healthy volunteers).
The presentation of work in this manuscript is separated into three main sections: simulation and associated optimization of the proposed technique, a comparison of the proposed method with existing methods via simulation, and an in vivo comparison in healthy volunteers.Fig. 1 gives a schematic overview of the workflow for reference.

Forward signal model
As per (Woolrich et al., 2006) the MR signal model used here describes the total signal (blood and static tissue).The model consists of an ASL component (which modulates the voxel's magnetization) and a BOLD component (contained within ΔR 2 ⁎ ) as described by Eq. (1).
where p n is the MR signal for n = 1 to N acquisitions, S n is the magnetization in a voxel, R 2 ⁎ , 0 is the baseline R 2 ⁎ , ΔR 2,n ⁎ is the change in R 2 ⁎ for each acquisition, and TE is the echo time.
As shown in Appendix A, p n can be estimated from a set of physiological parameters, end-tidal respiratory recordings, an ASL signal model, and a generalised BOLD calibration model.For a dual-echo ASL sequence, which is increasingly being acquired in calibrated fMRI experiments (Germuska and Bulte, 2014;Wise et al., 2013), the acquired data can be expressed as shown by Eq. ( 2).
where Y is the measured data, g is the multi-echo signal model, e is additive noise, θ is the set of physiological model parameters (see Table 1), and traces consist of P ET O 2 and ΔP ET CO 2 (measured in mm Hg).
Compared to previous implementations of dcfMRI, the calibration parameter is renamed (K), as it no longer includes contributions from OEF.Additionally, CVR is used as a replacement for fractional flow change using the empirical observation that ΔPaCO 2 is linearly related to ΔCBF for moderate changes in PaCO 2 (as caused by hypercapnic gas challenges) (Reivich, 1964;Tancredi and Hoge, 2013).
The forward model is solved using a regularized non-linear least squares fit to the acquired data, using a tuning parameter λ to determine the strength of the regularization (see Eq. ( 3).
The regularization is only applied for the first three parameters (K, OEF 0 and CVR), as these are the main drivers of any instability in the fitting.Because each parameter is centred and scaled prior to fitting, the relative strength of regularization is determined by our expectation of the biological variance of each parameter.To estimate the biological variance we have assumed a uniform distribution bounded by plausible in vivo values.Thus K (which is determined by the resting CBV) is assumed to lie between 0 and 0.3, OEF 0 is expected to lie between 0.1 and 0.7, and CVR is assumed to have a value between 1% and 6% ΔCBF/mm Hg.
Compared to previous analysis approaches the forward model estimates one additional parameter.Using a standard analytical approach to calculate CMRO 2 , estimates would need to be made of the BOLD signal changes due to O 2 and CO 2 , the calibration factor M, the change in blood flow due to CO 2 , and the resting blood flow.The extra parameter estimated in the forward model allows the CBF contribution to TE 2 (or any additional echoes) to be included into the analysis rather than averaged out (as is normal practice).The end-tidal traces (P ET O 2 and P ET CO 2 ) are used to constrain the fitting to the experimental paradigm.In addition, P ET O 2 is used to estimate the arterial change in blood T 1 as described in Appendix A.

Modelling studies
The proposed analysis method, and previous analysis methods (Bulte et al., 2012;Gauthier et al., 2012), were investigated with a noise modelling study.All data were generated using MATLAB, with the analysis being performed using a combination of MATLAB and FSL scripts.
Simulated data were generated using the described ASL signal model and a multi-compartmental BOLD model (Griffeth and Buxton, 2011).The multi-compartmental BOLD model extends that used for the analysis by explicitly considering contributions from extravascular and intravascular compartments for arterial, venous and capillary blood.Whereas, the BOLD model used for analysis groups all of these effects into three summary parameters (α, β and K).As such, the compartmental model is expected to be a more accurate description of the physiological dependences of the BOLD signal and therefore is utilised for data generation in this study.
Additive noise was modelled with an explicit physiological noise model (Kruger and Glover, 2001), which includes contributions from thermal white noise (σ 0 ), non-BOLD physiological noise (σ NB ), and BOLD-like physiological noise (σ B ).The Gaussian noise variances were modelled for 3.0 T as per Kruger et al., with σ NB = 0.21 and σ B =0.51 (TE = 29 ms) and σ B = 0.05 (TE = 2.7 ms).The voxel size was set to 3.4 × 3.4 × 7 mm, with a 6 mm 3D Gaussian smoothing kernel, producing σ 0 = 0.10.In addition, temporal correlations, using an autoregressive process, and low frequency drifts were included.The temporal AR(1) coefficient was set to 0.23 for the first echo and 0.55 for the second echo (Woolrich et al., 2006).Low frequency baseline drifts were estimated by fitting a 4th order polynomial to historic resting-state data acquired with a dual-echo PASL sequence with echo times of 2.7 ms and 29 ms (matching those used in the simulation).Because the acquired resting-state data was only 12 min in duration, interpolation in the frequency domain was used to artificially extend the data to the required time (18 min).
The simulated gas paradigm (see Fig. 2) was designed to match the end-tidal traces recorded in vivo.The simulated paradigm was 18 min in length with 3 periods of hypercapnia and 2 periods of hyperoxia interleaved around baseline/resting periods.For each hypercapnic period ΔPaCO 2 was set to 11 mm Hg and for each hyperoxic period ΔPaO 2 was set to be 240 mm Hg (values chosen to match in vivo recordings).
During periods of hypercapnia a 24-mm Hg increase was included in PaO 2 , while during hyperoxia a 2-mm Hg reduction in PaCO 2 was included.These interaction terms are based on the observed in vivo changes and it is assumed that the reduction in PaCO 2 during hyperoxia is responsible for any flow changes during these periods (i.e.there is no separate effect of hyperoxia on flow included in the modelling).
The simulated T 1 of arterial blood was set to match that used in the forward model (see Appendix A).A gamma-variate function (Madsen, 1992) was convolved with the generated block gas designs to account for the dispersion in the arterial stimulus.The t max and dispersion parameters of the gamma-variate were found from recordings of endtidal CO 2 and O 2 during in vivo scanning, and were set to t max =25s and dispersion = 0.5 for both the CO 2 and O 2 stimuli.While arterial changes in CO 2 are generally quicker than changes in O 2 ,inourimplementation O 2 stimuli are initiated with a period of 100% O 2 and followed up with a period of inspiratory hypoxia to ensure fast transitions.Thus, a similar transition time was found for both stimuli.
Simulations were undertaken to optimise the tuning parameter λ used in the forward modelling approach.One thousand time-courses were simulated with a set of randomly sampled physiological values for each simulation.The simulations were repeated with values of λ ranging from 0.25 to 5. The physiological values were sampled from a wide range to encapsulate potential in vivo variation both in healthy and diseased tissue (see Table 2 for a summary).
The simulated data was analysed with the forward model to find the optimal value of λ.The root mean squared error between the estimated and simulated OEF 0 for each value of λ was used to pick the optimal value.
The same simulated data was used to compare the performance of the forward model (with optimised λ) to the analytical methods proposed by Bulte et al. and Gauthier et al. (Bulte et al., 2012;Gauthier et al., 2012).For the analytical methods the simulated data was processed using FEAT (Jenkinson et al., 2012) to make estimates of ΔCBF| CO2 and ΔBOLD| CO2 and ΔBOLD| O2 .The change in blood flow and BOLD signals was analysed using a GLM fit to the data as per (Gauthier et al., 2012).However, rather than discarding the transition data (between baseline and the stimulus steady state) regressors were convolved with a gamma-variate function as per (Bulte et al., 2012).This was done to allow the best match between the forward model and the analytical approaches, both using all the time course data    from the CBF and BOLD acquisitions.For consistency with the forward model, oxygen and carbon dioxide are assumed to be isometabolic.However, a CBF reduction during O 2 was included separately from the model fits as per (Gauthier et al., 2012).A 2.4% reduction in flow was used, assuming a linear scaling of the flow reduction with the change in end-tidal O 2 .OEF estimates for the Bulte method were calculated by first calculating the CO 2 M value and then feeding this into the O 2 calibration equation, while the OEF for the Gauthier method was calculated using a bisection fitting method to find the values of M and OEF that fit the estimated CO 2 and O 2 signal changes (Bulte et al., 2012;Gauthier et al., 2012).

In vivo scanning (gas delivery)
In vivo data were acquired from a cohort of 10 healthy volunteers (6 male and 4 female, mean age 27.4 ± 10 years) with each subject having given informed written consent.Experimental procedures were approved by the local institutional ethical review committee.Respiratory challenges were administered using a previously described dynamic end-tidal forcing system (Wise et al., 2007).However, the feedback mechanism was disabled, producing an automated delivery system for fixed inspired fractions rather than targeting end-tidal values.Gas delivery was controlled by four mass flow controllers (MKS Instruments, Wilmington, MA, USA).Tidal gases, sampled from the volunteer's facemask and P ET CO 2 and P ET O 2 , were measured using rapidly responding gas analysers (AEI Technologies, Pittsburgh, PA, USA).The gas mixtures were delivered at a total flow rate of 25 l per minute (lpm) through a fast gas mixing chamber connected to a tightfitting facemask worn by the volunteer (Quadralite, Intersurgical, Wokingham, Berkshire, UK).The facemask was connected by a series of one-way valves to a reservoir limb (as described by Tancredi et al.   (2014)), which provided a gas reserve in case the instantaneous inspiratory rate exceeded 25 lpm.
The respiratory paradigm matched the design used in the simulation studies (Fig. 2).The design is 18 min in length and consists of 3 periods of hypercapnia interleaved with 2 periods of hyperoxia.Each period of hyperoxia was initiated with a short duration (14 s) of 100% oxygen, which acts as a pre-emphasis to shorten the transition from baseline to stable hyperoxia.At the end of each hyperoxic block a hypoxic gas mixture (10% O 2 , 90% N 2 ) is delivered for (40 s) to ensure a rapid transition back to baseline.Inspired fractions of 5% CO 2 and 50% O 2 with balanced air were delivered during plateau periods.

In vivo scanning (MRI acquisitions)
Image data were acquired on a 3 T whole body MRI system (GE Excite HDx, Milwaukee, WI) using an eight-channel receive-only head coil.Functional data were acquired using a pulsed arterial spin labelling (ASL) proximal inversion and control for off-resonance effects (PICORE), quantitative imaging of perfusion using a single subtraction (PICORE QUIPSS II) (Wong et al., 1998) imaging sequence.This sequence used a dual-echo gradient echo (GRE) readout and spiral k-space acquisition (Glover, 1999) (TE1 = 2.7 ms TE2 = 29 ms, TR = 2.2 s, flip angle 90°, FOV 22 cm, matrix 64 × 64, 12 slices of 7 mm thickness with an inter-slice gap of 1 mm, TI1 = 700 ms, TI2 = 1500 ms for the most proximal slice, 10 cm inversion slab thickness, adiabatic hyperbolic secant inversion pulse, 10 mm gap between labelling slab and bottom slice, 10 cm QUIPSS II saturation band thickness).Two scans were performed in each session with 490 repetitions (image volumes) collected in each scan.A separate single shot (M 0 ) scan was acquired with the same parameters as the functional runs to measure the equilibrium brain tissue magnetisation for purposes of calculating CBF 0 .AwholebrainT 1 -weighted structural scan (fast spoiled gradient recalled echo, 1 × 1 × 1 mm voxels, TI/TR/TE = 450/7.8/3ms) was acquired to facilitate segmentation of the grey matter and CSF.

Data analysis
To enable quantification of CBF 0 the signal from CSF was measured from the fully relaxed scan to estimate the M 0 of blood.Calculation of all physiological parameters was performed to match the pipeline used for analysing the simulated data (for both the forward modelling approach and the two analytical methods).Thus, spatial smoothing using a 6 mm FWHM Gaussian low-pass filter was applied to all data and a high-pass filter was applied to TE 2 data.When using the forward modelling approach TE 1 data was high pass filtered using surround subtraction (see Appendix A for more details).
Anatomical T 1 -weighted images were registered to the fMRI data using the FSL program FLIRT (Jenkinson and Smith, 2001).Partial volume estimates of grey matter (greater than 50%) (FAST image segmentation) were used to generate grey matter parameter estimates.The code used in this manuscript is openly available from the Cardiff University data archive http://dx.doi.org/10.17035/d.2015.100126.However, due to ethical considerations open access cannot be given to the in vivo subject data or data derived from this.
Statistical analysis of the results was performed using Matlab.Statistical tests include regression analysis, t-tests and Westenberg scale test (Westenberg, 1948) for interquartile range equality.In all cases a significance level of p = 0.05 was used as a threshold to determine statistical significance.

Modelling results
The noise modelling simulations demonstrated that a tuning parameter (λ) close to 1 produced the smallest RMS error in OEF 0 estimates Fig. 8.In vivo maps of OEF for a single subject (four slices).The top four slices are calculated using the forward modelling method; the bottom four slices are calculated using the method proposed by Gauthier et al.The maps generated with the forward model are more uniform with significantly fewer regions with non-physiological estimates.when using the forward model (see Fig. 3).Thus, this value was used in all subsequent modelling and data analysis.
The performance of the forward model was compared to analytical analysis of the data as per (Bulte et al., 2012;Gauthier et al., 2012)in both simulated and in vivo data.The simulated results demonstrate a significant reduction in uncertainty when using the forward model.Fig. 4 shows boxplots of the error in estimated simulated OEF 0 values for each of the three methods (the performance of the forward modelling approach is displayed with and without regularisation).The interquartile ranges (IQR) of the error in OEF 0 estimates were 0.15, 0.11, 0.18, and 0.17 for the forward model (without regularisation), the forward model (with L2 regularisation), the analytical approach proposed by Bulte et al. and the analytical approach proposed by Gauthier et al.The IQR of modelled OEF 0 estimates is significantly reduced when using the regularised forward modelling method (p b 0.05), as assessed by the Westenberg test for IQR equality.However, no significant difference was found between the two analytical approaches (p =0.06).
The median error in OEF 0 was − 0.009, − 0.010, 0.064 and 0.011 for the forward modelling approach, forward modelling with L2 regularisation, the Bulte analytical method and the Gauthier analytical method.No significant difference was found between the two forward modelling approaches (p = 0.785), suggesting that the inclusion of L2 regularisation has not introduced any bias into the analysis.However, a t-test revealed significant differences (p b 0.05) between each of the other methods, suggestive of a small bias inherent in each of the analytical methods.
In addition to the noise modelling results, estimates from the forward model were compared to the simulated data in a noise free simulation.These comparisons revealed that the resting venous blood volume could be estimated by a simple linear scaling of the calibration parameter K.This is consistent with the Davis BOLD model (as expressed by Eqs.(3)a nd(4)i nAppendix A), where K is defined as A ×CBV 0 .While the value of A will change with field strength and tissue structure, a typical value can be calculated by comparison with the multi-compartmental BOLD model.This assignment of a constant value to A is analogous to using fixed values for α and β, which is routinely done in calibrated BOLD analysis.Computing a linear regression of the fitted value of K against the modelled venous CBV gives a value of 3.7 for A with an R 2 of 0.78.

In vivo results
In all subjects end-tidal O 2 increased during periods of hyperoxia, while end-tidal CO 2 increased during hypercapnia.Fig. 5 shows an example end-tidal recording from one subject for O 2 and CO 2 .A s demonstrated by the traces, periods of hyperoxia appear to produce a reduction in P ET CO 2 , while periods of hypercapnia correspond with a slight increase in P ET O 2 .The average baseline P ET O 2 was 116 mm Hg and the average baseline P ET CO 2 was 43.5 mm Hg.Plateau levels of hyperoxia caused an increase of 240 mm Hg in P ET O 2 .The average change in P ET CO 2 from baseline was approximately 11 mm Hg.
Fig. 6 shows an example voxel-wise fit for TE 1 and TE 2 from an individual subject.Due to the high-pass filtering of TE 1 data any low frequency information is effectively removed from the time series (removing any baseline fluctuations and BOLD contributions).However, the presence of CBF information in TE 2 is evident.It thus highlights one of the advantages of the proposed method, in that it uses more of the available data to determine signal changes compared to a separate analysis of the BOLD and ASL data.BOLD and ASL signal changes in the TE 2 data are well matched to the fitted data, demonstrating the effectiveness of incorporating end-tidal traces in the fitting process.
The forward modelling results (with L2 regularisation) from each of the 10 individual subjects are shown in Table 3.As can be seen from the table, the median grey matter parameter estimates are w i t h i nt h ee x p e c t e dp h y s i o l o g i c a lr a n g ef o re a c ho ft h es u b j e c t s .The group mean OEF 0 estimate(acrossthesubjects)is0.39±0.03,with a mean venous CBV of 2.25% ± 0.36 and a mean CMRO 2 value of 200 ± 40 μmol/100 g/min.
Fig. 7 shows scatter plots of the relationship between the measured grey matter CBF 0 and the estimated grey matter OEF 0 for each subject calculated with each analysis method (the forward model and the two analytical methods).It is expected that the resting OEF should fall with increasing CBF 0 (Lu et al., 2008); this is clearly the case for the forward modelling approach with a strong negative correlation between the two parameters (R 2 =0.66,p b 0.05).However, no significant correlations are observed in either of the analytical approaches (p =0.46 and 0.57).The main reason for the lack of such a relationship appears to be the four subjects who exhibited a low resting CBF.These subjects have high OEF 0 estimates with the forward model, but low estimates with the analytical approaches.
Table 4 shows the estimated grey matter OEF 0 for each subject with each method, and the associated coefficient of variation (COV) for each estimate, where the COV represents the voxel-wise variation over space per subject.It is clear from the table that those subjects with the lowest resting CBF (1, 2, 7 and 10) have unusually low estimates of OEF 0 when using either of the analytical approaches.Paired t-tests show significant differences in OEF 0 estimates between each of the methods (p b 0.05) apart from the two analytical methods (p = 0.21).The coefficient of variation of OEF 0 estimates is significantly (p b 0.05) reduced when using the forward modelling approach, with an almost 2-fold reduction in the group average COV (0.35 compared to 0.65 and 0.67).A paired t-test also reveals that the COV for the Gauthier analytical method is larger than for the Bulte method.However, the size of this difference is very small.Fig. 8 shows example OEF 0 parameter maps calculated in an individual subject with each of the analysis methods (the forward model and the analytical approach).Consistent with the numerical assessment, the parameter maps of OEF 0 are significantly more uniform when calculated using the forward modelling approach.Additionally, parameter maps calculated with the analytical approach show several regions with areas of non-physiological estimates (OEF less than 0 or greater than 1).As shown by the numerical assessment, both of the analytical approaches produce nearly identical results; therefore we have only presented parameter maps from one method.The Gauthier implementation was chosen as the noise modelling study suggests this method is less prone to bias in OEF estimates.It should be noted that all parameter maps show grey matter voxels with a partial volume estimate greater than or equal to 30% (rather than the 50% threshold used in the numerical/statistical analysis).This cut-off has been chosen to provide a better visual assessment of parameter variation, and is especially valuable when visualising CBV parameter maps (Fig. 9), where estimates are expected to vary with grey matter content.
A summary image showing a single slice for each of the 10 subjects for each of the parameters estimated by the forward model is shown in Fig. 9.The maps demonstrate that OEF is for the most part uniform and unrelated to grey/white matter content.However, the maps of CMRO 2 , CBF and CBV show the expected grey/white matter contrast associated with blood volume variation.Maps of CVR are generally smooth although there are several regions of elevated (and suppressed) CVR.It is not clear if this is physiological in nature or an SNR artefact (CVR maps are expected to be the most affected by the limited SNR of the ASL data).

Discussion
We have demonstrated the measurement of OEF and absolute CMRO 2 using a forward signal model and compared the performance of the fitting method to two previously proposed analytical methods.The proposed method simultaneously estimates signal changes for TE 1 and TE 2 data (from a dual-echo acquisition) to find the physiological parameters that provide a best fit to the acquired data.Because all the model parameters are estimated simultaneously, information from CBF estimates is used to infer on BOLD signal changes and vice versa.Additionally, the whole time course data (for both echoes) is used in the analysis.This is different from previously proposed methods (Bulte et al., 2012;Gauthier et al., 2012;Wise et al., 2013), which either separate the acquired data into plateau periods or independently analyse the BOLD and CBF signal changes with a GLM (or use a combination of the two approaches).
In the proposed method a combined CBF-BOLD model is used in combination with end-tidal recordings to simultaneously estimate all BOLD and CBF changes from the acquired data.As with other numerical analysis techniques, the proposed method aims to make better use of the acquired data by considering all of the available information.In the implementation presented here we make use of regularisation to produce more stable estimates of the model parameters.As with any regularisation technique there is always a trade-off to be made between stability and accuracy of the results.In this work we have used a modelling study to determine the optimal value of the tuning parameter (λ)to avoid unnecessarily biasing estimates of OEF (and thus CMRO 2 ).The forward modelling method employed here enables a variety of numerical approaches to be used in tackling the minimization problem.Here we have presented a regularized least squares approach, which is a practical and computationally inexpensive method.However, such an approach is limited; for example, it makes implicit assumptions about the parameter and noise distributions.Alternatively, the problem could be considered in Bayesian framework where the regularization takes the form of prior distributions.Such an analysis method gives the opportunity to define more realistic noise models, or depending on inference method, reduce the strict distributional assumptions.If the reader is interested in pursuing a Bayesian approach the MATLAB code for a variational Bayes method is available http://dx.doi.org/10.17035/d.2015.100126.
The group mean OEF value found in this study, 0.39 ± 0.03, is similar to that found in previous dcfMRI experiments, which have reported values of 0.35 ± 0.04 (Gauthier et al., 2012), 0.36 ± 0.05 (Germuska and Bulte, 2014), 0.38 ± 0.14 (Bulte et al., 2012), and 0.40 ± 0.02 (Wise et al., 2013).However, previous dcfMRI analysis methods have reported a high number of voxels with invalid parameter estimates.For example, parameter maps presented in (Wise et al., 2013) reveal that estimates of OEF were unreliable in over half of the subjects for significant proportions of grey matter voxels.In addition, the results presented by Gauthier et al. (2012) depict numerous regions of invalid M and OEF estimates from individual subjects.The results presented by Bulte et al. (2012) do not explicitly address invalid parameter estimates.However, the reported individual OEF maps show large numbers of estimates close to zero.We have employed a combination of noise modelling and in vivo analysis to compare the stability of the proposed method with that of those previously proposed.The modelling results predict a significant advantage for the proposed method with regards to error in OEF estimates and bias.In addition, in vivo results demonstrate an increase in stability on a voxel-wise basis, with the proposed method being able to estimate OEF values in regions where previous methods were shown to fail.Furthermore, the relationship between the estimated in vivo parameters shows greater conformity to expected physiological variation with the proposed method.For example, data analysed with the forward modelling approach showed a strong negative correlation between resting CBF and OEF, which was absent when the data was analysed with previously proposed analysis methods.
Visual inspection of the estimated parameter maps also reveals the expected variation in CBV estimates with grey matter content.Such confirmation is important as it suggests correct functioning of the method at the voxel-wise level across a range of baseline physiology (in addition to the global correlations between CBF and OEF).Previous analysis methods have not always demonstrated this feature.For example, dcfMRI data analysed by Wise et al. (2013) reported maps of M that had a uniform appearance, rather than displaying the expected CBV modulation; while the individual OEF maps reported by Bulte et al. (2012), Gauthier et al. (2012) showed significant regional variation.
The modelling undertaken to compare the analysis methods was also employed to investigate the relationship between the simplified BOLD models routinely used in calibrated fMRI and the compartmental BOLD model used in the generation of synthetic data.It is apparent from the Davis BOLD model that a dcfMRI experiment is sufficient to determine the BOLD sensitive blood volume if the calibration factor A is known.By comparing the generated data (venous CBV) to the results produced with our forward model we found that A could be assigned a value of 3.7.In practice the method is similar to that proposed by Blockley et al., and produces similar estimates of venous CBV.The in vivo estimates of venous CBV in this work have a group mean grey matter value of 2.25 ± 0.36% and show the expected variation with grey matter content.These results are comparable to those reported in Blockley et al. (2013) 2.18%, An and Lin (2002) 2.46% and He and Yablonskiy (2007) 1.75%.
While the validity of the proposed method has been explored with modelling studies (and tested in vivo), one source of error, arising from the end-tidal traces, was not explicitly investigated.All of the BOLD models used in the analysis are developed from steady-state considerations and have not been fully explored for dynamic signal changes.In our design we have tried to minimise any contribution from this potential source of error by employing pre-emphasis to reduce the transition time between states.We expect such measures, as well as a tightly controlled experiment, to minimise these effects.However, care should be taken when employing these analysis methods when such fast transitions, or stable states, cannot be achieved.
In this paper we have developed a new analysis framework for dcfMRI.The proposed framework has been tested with noise modelling studies and in vivo acquisitions in healthy volunteers.The results from the in vivo analysis agree well with the expected physiology and demonstrate the improvements predicted by the modelling.The stability and reliability of the physiological estimates appear to be a significant improvement on existing techniques.We hope that these contributions will help move this emerging technique from a strictly research tool towards a clinically viable method capable of measuring resting OEF, CBV, CVR, and CMRO 2 from a single scanning session.simplified BOLD model (Wise et al., 2013).The ASL model is defined for a PASL QUIPSS II acquisition scheme (as is routinely used for ASL data acquisition at our institution).The BOLD model extends that proposed by Davis et al. (1998) and Hoge et al. (1999) to allow for simultaneous changes in both oxygen content and cerebral blood flow.This generalisation of the model allows the complete BOLD time series to be estimated for arbitrary levels of arterial oxygen content and alterations in CBF (provided isometabolic conditions can be reasonably assumed).The two models (ASL and BOLD) are combined, as described by Eq. ( 1), where the BOLD contribution is incorporated via a change in the transverse relaxation rate (ΔR 2 ⁎ ).
where p n is the MR signal for n = 1 to N acquisitions, S n is the magnetization in a voxel, R 2,0 ⁎ is the baseline R 2 ⁎ , Δ R 2,n ⁎ is the change in R 2 ⁎ for each acquisition, and TE is the echo time.
So that the BOLD signal model can be more easily incorporated into the ASL model it is expressed in terms of change in transverse relaxation rate (R 2 *), rather than the more common usage of change in signal magnitude, as shown below.

A.1.1. BOLD signal model
The contribution to the transverse relaxation rate due to deoxyhaemoglobin can be expressed as (Boxerman et al., 1995;Ogawa et al., 1993) where [dHb] is the deoxy-haemoglobin concentration, CBV is the BOLD sensitive blood volume, and A and β are scaling factors dependent on the vessel geometry and field strength.The change in relaxation rate at non-baseline values of CBV and deoxy-haemoglobin concentration can be expressed as where the subscript 0 represents the baseline value.
Rearranging Eq. ( 3) and exploiting the Grubb flow-volume relationship (Grubb et al., 1974) we get the following expression (which describes the change in transverse relaxation rate due to a change in blood flow and/or blood saturation): where α and β are assigned literature values (0.14 and 0.91 (Griffeth and Buxton, 2011) We are able to incorporate recordings of P ET CO 2 into Eq.( 4) by substituting the fractional change in blood flow for the relative cerebrovascular reactivity (CVR).This has two advantages; first it makes use of all the available information (incorporating the end-tidal recordings), second it allows parameter regularisation to be used that is not dependent on the experimental design (e.g.levels of CO 2 used as a stimulus).Thus Eq. ( 4) canbere-writtenas: where CVR is the flow response to CO 2 in %ΔCBF/mm Hg and ΔPaCO 2 is the change in end-tidal CO 2 from baseline in mm Hg.
The fractional change in deoxy-haemoglobin concentration can be expressed as follows (Wise et al., 2013), where again the change in CBF has been substituted for CVR.
where CaO 2 is the arterial oxygen content, which can be calculated from the end-tidal oxygen recordings.
As in previous dcfMRI studies (Bulte et al., 2012;Gauthier et al., 2012) the parameters α and β are assigned fixed values and assumed to remain constant throughout the acquisition period.Although K is not assigned a fixed value (as it represents as scaled version of the resting blood volume), it is assumed to remain constant during both hypercapnic and hyperoxic stimuli.

A.1.2. ASL signal model
As per Woolrich et al. (2006) the physiological ASL model includes a static compartment S n s and a blood compartment S n b such that the magnetization in a voxel is considered as a combination of the two, S n =S n s +S n b .The signal model for each tissue compartment is as described by Woolrich et al. unless otherwise noted.Thus, the magnetizations from the blood and static compartments are given by Eqs. ( 7) and ( 8) respectively: where M 0,blood is the resting arterial blood magnetisation per unit volume and R = Εe −TI2 /T1b .The tagging efficiency Ε is assumed to be 1.0 and δt is the transit delay (time between arterial tagging and the arrival of blood in the imaging region).The alternating tag and control conditions are represented by the term ρ n , which is 1 for control and −1fortag.
where M 0 is the baseline static magnetisation per voxel and ΔM is any CBF induced change in the magnetization.
Changes in M 0 are assumed to derive from a change in blood volume and VASO water exchange, as described by Eq. 33 in Woolrich et al. (2006).Thus, the hypercapnic change in S n s can be calculated from the resting blood flow (CBF 0 ), resting total blood volume (CBV t ) and the change in blood flow (using the flow-volume relationship described by Grubb et al. (1974).When simulating data all of these variables are known, and thus the change in M 0 can be simulated for a range of physiological conditions.However, CBV t is an unknown parameter when analysing data.Thus, we chose a pragmatic approach of using a fixed CBV t of 5% when analysing data, with the change in M 0 being estimated from the flow parameters.
In an extension to the model proposed by Woolrich et al., the T 1 of blood is allowed to vary throughout the time series.This is necessary in the context of dcfMRI to account for the T 1 effects of dissolved arterial oxygen during periods of hyperoxia.The relationship between PaO 2 and arterial blood T 1 is taken to be linear and fitted to the data presented in (Ma et al., 2014), as described by Eq. (9): where b is expressed in mm Hg, and the intercept c can be related to the normoxic T 1 of blood as follows c =1.725− b •110.The normoxic value at 3.0 T (1.725 s) is taken from the literature (Lu et al., 2004).
A fixed value of − 5×10 −4 was used for b, extrapolating from Ma et al. (2014), and thus a value of 1.78 for c.These values were used to convert PaO 2 to T 1,b for all voxel-wise analysis (both simulated and in vivo).

A.2. High-pass filtering
The effect of slow static tissue signal changes was minimised for TE 1 and TE 2 by fitting to high-pass filtered versions of the signals.For TE 1 high-pass filtering was carried out using surround subtraction, where y[n] TE1 is the acquired TE 1 data, h[n]=[−12 -1] / 2 and ŷ [n] TE1 is the high-pass filtered signal.The high-pass filter for TE 2 used a lower cut-off frequency (to preserve the BOLD information) and was implemented to match that provided by fslmaths (Jenkinson et al., 2012).The mean signal for each echo time is preserved so that R 2,0 ⁎ and M 0 can still be estimated.To minimise any bias introduced by the high-pass filters the forward modelled data was processed with matching filters.
Although changes in CBF during the time series will alter the transit delay (δt), this parameter is assigned a constant value (0.4 s).The use of a constant transit delay in the context of QUIPSS II ASL is considered reasonable as this acquisition scheme is designed to be insensitive to changes in δt (Wong et al., 1998).

Fig. 1 .
Fig. 1.Schematic overview of the key pieces of work presented in this manuscript.

Fig. 5 .
Fig.5.Example end-tidal gas recordings from an individual subject.Top: End-tidal oxygen recordings.Bottom: End-tidal carbon dioxide recordings.The rise and fall times for both stimuli appear to be similar.Periods of elevated P ET CO 2 produce a corresponding increase in P ET O 2 , while periods of elevated P ET O 2 appear to produce a reduction in P ET CO 2 .

Fig. 4 .
Fig. 4. Boxplot showing the error between the simulated and estimated OEF 0 for each analysis method. A. Forward model.B. Forward model with regularisation.C. Bulte analytical method.D. Gauthier analytical method.The forward modelling method with regularisation demonstrates the least bias and smallest error.The Bulte and Gauthier analytical methods demonstrate similar levels of uncertainty.There is a significant difference (p b 0.05) in the spread of OEF 0 estimates between each method, apart from the two analytical approaches (p = 0.06).A Student's t-test of the means shows no significant difference between the forward model with or without regularisation (p =0.785).However, each of the other methods has significantly different means (p b 0.05).

Fig. 3 .
Fig.3.RMS error in estimated OEF 0 vs. simulated OEF 0 for the forward modelling analysis method plotted against the value of the regularisation tuning parameter λ.The optimal value of λ (minimum error) is close to 1.

Fig. 7 .
Fig. 7. Scatter plots of resting CBF against resting OEF as estimated by each of the implemented analysis methods.There is a strong negative correlation between the resting CBF and OEF when the in vivo data is analysed with the forward model.No such relationship is found when analysing the data with the analytical approaches of Bulte and Gauthier.

Fig. 6 .
Fig. 6.Example voxel-wise fit of in vivo data.Top: Processed TE 1 data (blue) and fit(black).Bottom: TE 2 data (blue) and fit (black).Both the acquired and estimated TE 2 data have been high-pass filtered (300 s cut-off) to reduce the influence of baseline drift.

Table 1
Summary of the physiological parameters used in the forward signal model.

Table 2
Parameter values and ranges used in the simulation studies.

Table 3
Average (median) grey matter parameter estimates for each subject assessed using the proposed forward modelling analysis method.

Table 4
Average (median) grey matter parameter estimates and the coefficient of variation (COV) of OEF 0 estimates for each subject and analysis method (forward modelling with L2 regularisation and the two analytical approaches proposed byBulte and Gauthier).