Calibrated BOLD using direct measurement of changes in venous oxygenation

Calibration of the BOLD signal is potentially of great value in providing a closer measure of the underlying changes in brain function related to neuronal activity than the BOLD signal alone, but current approaches rely on an assumed relationship between cerebral blood volume (CBV) and cerebral blood flow (CBF). This is poorly characterised in humans and does not reflect the predominantly venous nature of BOLD contrast, whilst this relationship may vary across brain regions and depend on the structure of the local vascular bed. This work demonstrates a new approach to BOLD calibration which does not require an assumption about the relationship between cerebral blood volume and cerebral blood flow. This method involves repeating the same stimulus both at normoxia and hyperoxia, using hyperoxic BOLD contrast to estimate the relative changes in venous blood oxygenation and venous CBV. To do this the effect of hyperoxia on venous blood oxygenation has to be calculated, which requires an estimate of basal oxygen extraction fraction, and this can be estimated from the phase as an alternative to using a literature estimate. Additional measurement of the relative change in CBF, combined with the blood oxygenation change can be used to calculate the relative change in CMRO2 due to the stimulus. CMRO2 changes of 18 ± 8% in response to a motor task were measured without requiring the assumption of a CBV/CBF coupling relationship, and are in agreement with previous approaches.


Introduction
Blood oxygenation level dependent (BOLD) contrast is widely used in functional magnetic resonance imaging (fMRI) to monitor brain function. However, the BOLD signal does not provide a direct measurement of brain function, but rather it monitors the haemodynamic response to changes in the underlying brain function. As a result the BOLD signal is blurred both spatially and temporally with respect to the underlying changes in neuronal function, and the amplitude of the BOLD signal has a complex, indirect relationship on the amplitude of the neuronal activity, which depends amongst other things on the local vascular structure and reactivity (Boynton et al., 1996;Buxton et al., 1998;Friston et al., 1998;Turner, 2002;Shmuel et al., 2007). These vascular confounds restrict the conclusions that can be drawn from a BOLD experiment, particularly in situations where either the baseline cerebral blood flow (CBF) and cerebral blood volume (CBV), or the capacity of the cerebral vasculature to respond may be altered, such as in pharmacological studies or some pathologies. Baseline and dynamic haemodynamic properties will combine, along with the cerebral metabolic rate of oxygen consumption (CMRO 2 ), to determine the amplitude of a BOLD signal response to changes in brain function.
BOLD contrast depends primarily on changes in venous blood oxygenation (Y) and venous cerebral blood volume (vCBV) (Ogawa et al., 1993), and in turn, Y depends on CMRO 2 and CBF, whilst the vCBV is thought to be coupled to CBF (Grubb et al., 1974;Buxton et al., 1998). Thus the BOLD signal depends on the changes in vCBV, CBF and CMRO 2 .
CMRO 2 is closely related to tissue energy demand (Rothman et al., 1999), so is expected to provide a more direct measure of underlying brain function than the BOLD signal, as it is less affected by local haemodynamic properties. Task related changes in CMRO 2 have previously been calculated by measuring BOLD and CBF changes in response to a task (Davis et al., 1998). This approach requires knowledge of a calibration parameter, which is the BOLD signal change that would be achieved if the venous blood volume were fully oxygenated. This calibration parameter can be calculated using measured BOLD and CBF changes in response to hypercapnia (Davis et al., 1998), or by using BOLD changes in response to hyperoxia combined with a model estimating the effect of hyperoxia on venous blood oxygenation .
Regardless of whether hypercapnia or hyperoxia is used, this technique assumes that CBV is coupled to CBF for steady state conditions, instead of directly measuring CBV. This coupling relationship was initially measured in an animal model (Grubb et al., 1974) as a power law relationship, such that where α = 0.38. Results subsequently obtained in humans suggested that α = 0.29 (Ito et al., 2003). However, this relationship refers to the total CBV, whereas the BOLD signal arises primarily from the vascular component containing deoxygenated haemoglobin (dHb) (i.e. vCBV corresponding to veins and venules). Furthermore, vCBV has been shown to increase relatively less than total CBV, in response to both hypercapnia (Lee et al., 2001) and forepaw stimulation (Kim et al., 2007) in rats. The coupling relationship between vCBV and CBF has been measured in humans as α = 0.23 during a combined motor and visual task (Chen and Pike, 2009) and α = 0.18 in response to hypercapnia (Chen and Pike, 2010), indicating that using the total CBV rather than vCBV coupling relationship for calibrated BOLD will lead to an underestimation of CMRO 2 (Chen and Pike, 2009). Nonetheless the coupling between vCBV and CBF remains poorly characterised in humans, and may vary with brain region and in pathology. In this work we propose a new approach to BOLD calibration, using BOLD contrast on hyperoxia and normoxia, to calculate the relative change in venous CBV and Y due to the stimulus, avoiding the need to make any assumption about the coupling between vCBV and CBF.

Theory
Hyperoxia causes an increase in cerebral venous blood oxygenation (Yv), providing a BOLD contrast which has previously been used to measure both absolute resting vCBV  and the fractional change in vCBV in response to a task (Blockley et al., 2012). It will be shown that if the same task is performed at both normoxia and hyperoxia then the relative change in Yv due to the task can also be measured, assuming that hyperoxia has no effect on either CBF or CMRO 2 . If the change in CBF due to the task is also measured, then the relative change in CMRO 2 (rCMRO 2 ) due to the task can be calculated (Davis et al., 1998) without assuming the relationship between vCBV and CBF. Table 1 lists the parameters used in the following description of the method.

Model for calculating relative task-related CMRO 2 changes
An analytical model is used to describe the effect of hyperoxia on tissue R 2 *, (Yablonskiy and Haacke, 1994).
where R 2,0 * is the transverse relaxation rate of tissue containing only fully oxygenated blood vessels, Y is the blood haemoglobin oxygen saturation fraction and V is the volume fraction occupied by blood vessels. The term k is a constant arising from an extravascular signal model, based on the static dephasing regime of spins located around randomly orientated blood vessels: where the 4π/3 term describes the effect of vessels with random orientation, Δχ is the susceptibility of deoxygenated haemoglobin relative to tissue, [Hb tot ] is the total haemoglobin concentration and B 0 is the static magnetic field. This model assumes a linear relationship between R 2 * and deoxygenated haemoglobin fraction (β = 1 in Davis model (Davis et al., 1998)) which is more appropriate for 7 T (Yablonskiy and Haacke, 1994;Driver et al., 2010;Croal et al., 2012a). The venous dHb fraction Q is defined such that Q = (1− Yv) and Q 0 denotes the normoxic rest condition (which is also equivalent to oxygen extraction fraction (OEF)). Hyperoxia increases venous oxygenation Yv, such that on hyperoxia Q = (Q 0 + ΔQ h )= Q 0 (1 + q h ), where ΔQ h is the absolute change in venous dHb fraction due to hyperoxia and q h is the relative change in venous dHb due to hyperoxia. Note that since hyperoxia will cause an increase in Yv, ΔQ h and q h will be negative. An analogous term q act can be defined for the relative change in Q due to task related changes in CMRO 2 and CBF. Therefore, the venous dHb fraction during a task performed at hyperoxia is Q =(Q 0 + ΔQ h + ΔQ act ) =Q 0 (1 + q h + q act ), assuming that the amount of oxygen extracted by the capillary bed is independent of hyperoxia at steady state (i.e. q act is the same at normoxia and hyperoxia).
The transverse relaxation rates during the neuronal activation R 2,act * and rest R 2,rest * at a given level of hyperoxia (q h ), can be modelled as where V 0 is the resting vCBV and ΔV act is the absolute change in venous CBV due to neuronal activation. R 2,0 * is the transverse relaxation rate of the tissue neglecting any venous effects, and includes effects of shimming. The method for estimating q h is discussed in Estimating q h theory section. If R 2,act * and R 2,rest * are measured for two or more values of q h then the gradient of a plot of R 2,rest * versus (1+q h ) ( Fig. 1a) will be kV 0 Q 0 , which is the calibration constant M used in previous calibrated BOLD experiments (Davis et al., 1998;Chiarelli et al., 2007a). The gradient of the plot of R 2,act * versus (1+q h ) (Fig. 1a) is k(V 0 + ΔV act )Q 0 , which will be referred to as M′. The difference in intercepts of the rest and activation curves is k(V 0 + ΔV act )Q 0 q act , which can Absolute change in Q due to the motor task q act Relative change in Q due to the motor task (such that ΔQ act = Q 0 ·q act ) CBF Cerebral blood flow rCBF Relative change in CBF due to the motor task (rCBF = ΔCBF/CBF 0 ) CBV Cerebral blood volume vCBV Venous cerebral blood volume rvCBV Relative change in vCBV due to the motor task (rvCBV = ΔvCBV/vCBV 0 ) CMRO 2 Cerebral metabolic rate of oxygen consumption rCMRO 2 Relative change in CMRO 2 due to the motor task (rCMRO 2 = ΔCMRO 2 /CMRO 2,0 ) α Coupling relationship between CBV and CBF (see Eq. (1)) β Power law relationship relating extravascular transverse relaxation to Y M Calibration parameter, as used in previous BOLD calibration studies be divided by M′ to give q act . The ratio of the activation to rest gradients gives the relative venous CBV (rvCBV) change on activation: If the relative change in CBF (rCBF) on activation can be measured (e.g. using arterial spin labelling (ASL)), then Fick's principle can be combined with the value of q act to estimate the relative change in CMRO 2 (rCMRO 2 ): To summarise, this method does not assume a coupling relationship between CBV and CBF (Grubb et al., 1974), but does assume negligible vasoconstriction due to hyperoxia (which is reasonable in the case of normocapnic hyperoxia (Croal et al., 2012b)) and that oxygen consumption is independent of hyperoxia so that q act is the same at normoxia and hyperoxia (an area that needs further investigation). This approach also assumes that the haematocrit remains constant both during hyperoxia and the motor task. Furthermore to create the plots described above, q h , the relative change in venous dHb due to hyperoxia, must be known; a method for determining q h is described in the next section.

Estimating q h
During hyperoxic periods, arterial oxygen partial pressure (PaO 2 ) is increased, with most of this extra oxygen being dissolved in arterial plasma, since arterial haemoglobin is already close to being fully saturated at normoxia. Oxygen dissociation in blood is described by a widely accepted model (Severinghaus, 1979). At 310 K and a pH of 7.4, the relationship between arterial haemoglobin oxygen saturation (SaO 2 ) and PaO 2 is plotted in Fig. 1b. Here end-tidal PO 2 (P ET O 2 ) was used to estimate PaO 2 (i.e. P ET O 2 = PaO 2 ) assuming that arterial blood gases are in equilibrium with alveolar gas. This assumption is not true in general as there are partial pressure gradients for oxygen between exhaled gas and alveolar gas due to non-uniform alveolar distribution of inhaled gas, and between alveolar gas and arterial blood due to the non-uniform distribution of lung perfusion (Ayres et al., 1964). The sum of these effects results in end-tidal to arterial gradients of less than 50 mm Hg in healthy subjects (Ayres et al., 1964). This PO 2 difference at hyperoxia would correspond to an error in calculation of plasma O 2 content of 0.15 ml per decilitre, or less than 0.7% of total blood O 2 content which would have a negligible effect on the calculation of q h . However in our study, end-tidal to arterial PO 2 gradient was reduced further by several factors unique to our method of controlling the PaO 2 (see Respiratory paradigm section). The employment of the sequential gas delivery circuit to administer the O 2 (in which gas in equilibrium with alveolar blood gas occupies the anatomical dead space) resulted in a more homogeneous distribution in alveolar PO 2 (Swenson et al., 1994;Brogan et al., 2004;Ito et al., 2008). This minimises the gradient between expired and alveolar gas, and also eliminates the effect of the variability of distribution of pulmonary blood flow to PaO 2 (Ito et al., 2008;Fierstra et al., 2011). The total arterial blood oxygen content (CaO 2 ) is the sum of the oxygen bound to haemoglobin and that dissolved in plasma where ϕ = 1.34 ml(O 2 )/g is the oxygen carrying capacity of haemoglobin, [Hb] = 15 g/dl blood is the concentration of haemoglobin and ε = 0.0031 ml/(dl blood ·mm Hg) is the solubility coefficient of oxygen in blood . Dissolved oxygen meets part of the tissue oxygen requirement (Rostrup et al., 1995). Venous haemoglobin oxygen saturation (SvO 2 ) can be estimated using the method proposed by Chiarelli et al. (2007a). The total venous blood oxygen content (CvO 2 ) is what remains after oxygen has been extracted from capillaries, and is given by where OE is the oxygen extraction. OE is assumed to be independent of the level of hyperoxia, and OE = CaO 2,0 ⋅ OEF, where CaO 2,0 is the total arterial oxygen content at normoxia and OEF is the oxygen extraction fraction at normoxia (= Q 0 ). Since venous haemoglobin is not close to being fully saturated under normobaric hyperoxia, it can be assumed that a negligible amount of the CvO 2 is carried as dissolved oxygen. This means that Eqs. (9) and (10) can be combined to give venous oxygen saturation: which is also plotted in Fig. 1b. This model can be used to estimate (1 − Q) from PaO 2 , which can be estimated by measuring exhaled gas oxygen fraction. Since OE is assumed to be independent of hyperoxia level, ΔQ h can be estimated independently of OE or OEF, but Q 0 must be estimated. There are two possible ways of doing this and we have attempted both in this paper. First the value can be taken from literature values of OEF, (we used 0.4 here). Alternatively, in principle it can be measured from venous blood (e.g. from the T 2 (Lu and Ge, 2008) or susceptibility as used here (Jain et al., 2010), see Analysis section).

Materials and methods
Nine healthy volunteers (3 male, 6 female; mean age = 27 ± 3 years, range 23-30 years) participated in this study. Ethical approval was given by the University of Nottingham Medical School Ethics Committee and all subjects gave informed written consent prior to participating. For one subject there was a significant change in P ET CO 2 (>1 mm Hg) during the transitions into both hyperoxia periods, which overlapped with the motor trials at hyperoxia and so this subject was discarded from further analysis.

Motor paradigm
All subjects were asked to perform a motor task during a respiratory challenge, so that the same motor task was repeated both at normoxia and hyperoxia. The motor task consisted of a bilateral sequential finger tap. This was visually cued, with the word 'TAP' displayed in red on a dark background on a projector screen at the end of the magnet bore. This was replaced by a white '+' fixation point during rest periods.
The respiratory challenge consisted of two repeats of 3 min of normoxia, followed by 3 min of hyperoxia, with a final 3 min of normoxia at the end. The motor paradigm consisted of blocks of two trials of 30 s ON/30 s OFF at each gas level, with the block finishing just before the start of a transition between the two gas levels (see Fig. 2A -Paradigm A). For the final subject, the paradigm was modified to maintain the subject's attention better: evenly spaced motor trials of 90 s of rest followed by 30 s ON/60 s OFF were superimposed on 2 min of normoxia, followed by two repeats of 2 min hyperoxia and 4 min of normoxia (see Fig. 2B -Paradigm B). Alternate trials fell at a respiratory transition and could not be used, giving 4 usable trials in total.
The last 4 subjects (indicated in Table 3) were also asked to complete an additional motor task during an ASL acquisition to allow the relative change in CBF to be measured, and thus the relative change in CMRO 2 to be estimated. This task consisted of 10 trials of Paradigm A or 8 trials of Paradigm B, performed with the subject breathing medical air.

Respiratory paradigm
A feed-forward, low gas flow system (RespirAct™, Thornhill Research Inc., Toronto, Canada) and a sequential gas delivery (SGD) breathing circuit (Banzett et al., 2000;Slessarev et al., 2007) were used to target end-tidal PCO 2 (P ET CO 2 ) and PO 2 (P ET O 2 ) independently (Slessarev et al., 2007). Source gases used by the system were O 2 , air, and two gas blends of N 2 , CO 2 and O 2 , so that all source gases contained safe concentrations of O 2 . The RespirAct™ follows the approach of Slessarev et al. (2007) to calculate the required flows of these source gases into the SGD breathing circuit to attain the targeted P ET CO 2 and P ET O 2 values. For the normoxic condition, both P ET O 2 and P ET CO 2 were maintained at the subject's resting values (~110 mm Hg and~40 mm Hg, respectively). For the hyperoxic condition, P ET O 2 was targeted at 500 mm Hg, whilst P ET CO 2 was maintained at the resting value.

Data acquisition
Scanning was performed on a Philips Achieva 7 T system, with head volume transmit and 16 channel SENSE head receive coil. Gradient echo (GE) EPI data were acquired every 2.4 s throughout the respiratory and motor tasks. Images consisted of 2 mm isotropic voxels, with a 192 × 192 mm 2 field-of-view, and 20 axial slices (2 mm thickness, no slice gap) spanning the motor cortex. Imaging parameters were TE = 25 ms, SENSE factor = 3, voxel band-width = 41.5 Hz, TR = 2.4 s and flip angle = 75°. ASL datasets were acquired using a FAIR labelling scheme, with TI = 1400 ms, selective thickness 10 mm wider than the imaging volume, non-selective thickness of 300 mm and background suppression pulses at 402 and 639 ms (Garcia et al., 2005); in-plane pre-and post-saturation were used. These images consisted of 2 × 2 × 4 mm 3 voxels, with a 192 × 192 mm 2 field-of-view, and 8 axial slices (4 mm thickness, no slice gap) with the same orientation and centre as the GE EPI data. Imaging parameters were TE = 14 ms, SENSE factor = 3, voxel band-width = 41.5 Hz, TR = 3 s (6 s for a tag/control pair). Two equilibrium magnetisation images were acquired for signal normalisation (same parameters as for the ASL, except with no inversion and a long TR of 10 s), one before and one after the ASL acquisition.
Analysis GE EPI datasets were motion corrected using MCFLIRT (FSL, fMRIB, Oxford, UK). Voxelwise linear detrending was then performed using a linear fit (MATLAB, The MathWorks, Natick, USA) to baseline timepoints (i.e. during periods of normoxia and not within 30 s of a previous finger tap). Linear detrending was chosen in preference to a high-pass temporal filter, due to the long cycle lengths of the paradigm. Maps of statistical significance of the BOLD response to the motor task were formed using the FEAT (FSL, fMRIB, Oxford, UK) general linear model. This was done on the datasets before voxelwise normalisation, but after motion correction and linear detrending. A boxcar design, representing the hyperoxia response was included in the model design as an independent regressor, so the hyperoxiabased BOLD response could be separated from the motor response in the analysis. 'BOLD motor activation masks' were formed of clusters of voxels with Z > 5 and P cluster b 0.05 for the response to the motor task.
Since the ASL data was background suppressed, this made motion correction problematic due to the low signal of the tag and control images. Instead the two equilibrium magnetisation images from the start and end of the ASL data set were subtracted to determine whether there was any systematic motion through the dataset, and no more than 1 voxel of displacement was found across all subjects. The interleaved acquisition of tag and control ASL images meant that linear temporal interpolation (MATLAB, The MathWorks, Natick, USA) had to be performed on both the tag and control datasets. ASL subtraction was then performed between interpolated tag and control datasets, resulting in a CBF-weighted timecourse with timepoints every 3 s, but a real temporal resolution of 6 s. Maps of statistical significance of the CBF response to the motor task were formed using the FEAT (FSL, fMRIB, Oxford, UK) general linear model. In this analysis, a high-pass filter was used with a cut-off period of twice the trial length, corresponding to 120 s for Paradigm A and 180 s for Paradigm B. 'CBF motor activation masks' were formed of clusters of voxels with Z > 2.3 and P cluster b 0.05. The lower Z-statistic threshold was used due to the intrinsically lower SNR of ASL compared with GE EPI BOLD.
Timecourses of measured P ET O 2 were used to estimate the absolute change in venous deoxyhaemoglobin fraction due to hyperoxia (ΔQ h ) (see Theory section). Q 0 and hence the relative change (q h ) was then calculated using an assumed resting OEF of 0.4. In addition, in an attempt to avoid the need to assume resting OEF, phase data constructed from the GE EPI dataset were used to calculate an alternative value of Q 0 and hence q h . The phase data were unwrapped in space over each volume (PRELUDE, FSL) and then unwrapped over time (UNWRAP, MATLAB). The motion correction transformations from the magnitude data were applied to the unwrapped phase data (MCFLIRT, FSL). To remove the global phase shift caused by increased oxygen concentrations in the frontal sinus and nasal cavity (Driver et al., 2011;Pilkinton et al., 2011), the phase data were high-passed filtered using a homodyne filter with a smoothing kernel with a FWHM of 4 mm (Hammond et al., 2008). The filtered phase data were then separated into normoxia and hyperoxia periods, defined as the final 2 min of each period (1 min for paradigm B). A line profile was taken across the sagittal sinus, to estimate the extravascular phase shift due to dHb in the sagittal sinus both at normoxia and hyperoxia (example in Fig. 3). A ratio of hyperoxia: normoxia phase (a) was calculated by a least-squares fit of the line profiles to the equation ϕ hyperoxia = a·ϕ normoxia . Voxels with a large intravascular component (such as sagittal sinus) were excluded from the fit, by applying a high intensity threshold to the magnitude data, averaged over the normoxia period. Assuming that the susceptibilities of tissue and blood plasma are equal to that of water and that a is also equal to the ratio of hyperoxia:normoxia susceptibilities, relative to tissue; the following equation can be used to relate a and ΔQ h to Q 0 (see Appendix A): where Δχ oxy = − 0.017 × 10 −6 , Δχ deoxy = + 0.247 × 10 −6 (Spees et al., 2001) are the volume susceptibilities of oxygenated haemoglobin and deoxygenated haemoglobin relative to water, respectively. For the GE EPI datasets each voxel's timecourse was normalised by dividing by the average value over all baseline timepoints and then average timecourses were calculated over all voxels within the 'BOLD motor activation mask'. In some cases, transient changes in P ET CO 2 occurred during a P ET O 2 transition, which persisted into the start of the next motor trial. Any motor trial which included a change in P ET CO 2 of greater than 1 mm Hg compared to the normocapnia value was discarded. The average signal during rest and active periods (motor task) were compared to the normoxic baseline (defined earlier), to estimate rest and an active % signal change (%BOLD) for each trial and each hyperoxia/normoxia condition. These % changes were converted to ΔR 2 * changes using the approximation ΔR 2 * ≈− %BOLD/(100TE), to estimate q act . Linear regressions (MATLAB, The MathWorks, Natick, USA) were performed to fit ΔR 2 * as a linear function of q h , separately for the rest and active conditions. The fitted parameters were then used to estimate q act and rvCBV, as proposed in the Theory section.
For each of the four subjects for which ASL data were acquired an additional 'combined mask' was formed from the intersection of the 'BOLD motor activation mask' and 'CBF motor activation mask'. The CBF-weighted data (after interpolation and subtraction) were normalised to the rest periods, by dividing by the average of the last half of all rest periods, and an average timecourse was calculated over all voxels in the 'combined mask'. The % change in CBF on activation was calculated by averaging over the active (motor task) period of all trials, multiplied by 100. Values of q act were recalculated for the 'combined mask' using the same method as described above, rCBF and q act were used to calculate rCMRO 2 , using Eq. (7). Maps of the pixel values of M, rvCBV and q act were also calculated.

Results
Transitions in ΔP ET O 2 of 330 ± 20 mm Hg (mean ± SEM over subjects) were achieved, which corresponded to an increase of ΔYv = 0.068 ± 0.003 and an 11% increase over the resting normoxia value. The average change in end tidal CO 2 on hyperoxia compared to normoxia was ΔP ET CO 2 = − 0.5 ± 0.1 mm Hg (range of − 1.1-0 mm Hg across subjects) which was assumed to be too small to have a significant effect on CBF and the resulting BOLD signal. Fig. 4 shows example BOLD, P ET O 2 and P ET CO 2 timecourses for subject 3 in the 'BOLD motor activation mask'. Averaged across all subjects, the baseline BOLD signal increased by 6.4 ± 0.9% with hyperoxia and the BOLD response to the motor task was increased by 37 ± 7% on hyperoxia compared to normoxia. An example of an average trial of the motor task at normoxia and hyperoxia, as well as an example of the linear fits are shown in Fig. 5. The gradients of the linear fits gave M = 36 ± 5% and M′ = 48 ± 7%. Individual subject values for the relative increase in vCBV, relative decrease in deoxyhaemoglobin fraction q act and corresponding absolute increase in ΔY act (= − Q 0 · q act ) during the motor task are reported in Table 2. Monte Carlo simulations were performed to estimate the precision of the rvCBV and q act values; 10,000 timecourses were simulated with the same temporal SNR (161) as the actual data and the resulting absolute standard deviation in the estimated values was σ(rvCBV) = 4.3 % and σ(q act ) = 2.1 % with no significant systematic error.
An example CBF average trial timecourse is shown in Fig. 5c. Results based on the combined mask, rCBF, q act and rCMRO 2 are reported in Table 3, for both an assumed Q 0 value of 0.4 and for Q 0 calculated from the phase data. For the 'combined mask', values of M and M′ were calculated as M = 28 ± 2% and M′ = 34 ± 2%, and rvCBV = 22 ± 7% was lower in all four subjects compared with those values calculated from the 'BOLD motor activation mask'.
Voxel-by-voxel maps of M, rvCBV and q act were formed by performing the above analysis on each voxel in the 'combined mask' and are shown in Fig. 6.

Discussion
This paper has described a new method of using hyperoxia to measure CMRO 2 . By comparing the response to a task at hyperoxia and normoxia it provides an extra degree of information on the vascular nature of the BOLD signal. As a result, this method does not make any assumption about the coupling between CBV and CBF, but instead provides a direct measurement of the change in blood volume and oxygenation on activation with relatively low sensitivity to noise as indicated by the Monte Carlo simulations. Essentially this method is equivalent to the previous method proposed for hyperoxia based measurement of CMRO 2  except that instead of using the Grubb relationship, here rvCBV is measured directly (Blockley et al., 2012). Recently, hypercapnia and hyperoxia calibrated BOLD approaches have been combined to map basal OEF (Bulte et al., 2012;Gauthier and Hoge, 2012), relying on an assumed coupling relationship between CBV and CBF. These basal measurements could be extended by including a hypercapnia challenge repeated both at normoxia and hyperoxia to overcome the need to assume this relationship.
The rCMRO 2 changes measured here agree with measurements made using hypercapnia-based calibration for a bilateral finger tap, with rCMRO 2 =16 ±9 % (Kastrup et al., 2002) and 5-45 % , and for a right-handed finger tap using positron emission tomography giving rCMRO 2 = 11 ±13 % (Ito et al., 2005). The value of rvCBV = 22 ± 7% measured here is approximately twice that measured for a motor task by Chen and Pike (2009). Possible reasons for this difference could be due to differences in task performance, the region selected or the methods themselves, and it would be interesting to compare the two methods directly in future work. In this study there was a trend for a smaller increase in rvCBV in the 'combined mask' defined by both BOLD and CBF activated regions, compared with the 'BOLD mask'. This suggests a smaller change in rvCBV in the microvasculature, compared with draining veins that contribute to the 'BOLD mask' but probably do not represent active tissue (Leontiev et al., 2007). The 'combined mask' was used for the estimate of CMRO 2 , to help to meet the extravascular signal assumption, by excluding large draining veins, whilst focusing the ROI on areas of active tissue. This paper reports the first measurements of M in humans at 7 T. Equivalent values measured in the motor cortex are M ≈ 4 − 6 % at 3 T using hyperoxia or hypercapnia Chiarelli et al., 2007b;Mark et al., 2011) and M ≈ 5 − 9 % at 1.5 T using hypercapnia (Kastrup et al., 2002;Stefanovic et al., 2004Stefanovic et al., , 2005Stefanovic et al., , 2006). An alternative approach at 3 T, calculating M by attempting to fully saturate venous blood using combined hypercapnia and hyperoxia (10% CO 2 carbogen) reported M = 7.5-9.5% in the visual cortex (Gauthier et al., 2011). The increased value reported here (M = 36 ± 5 %, 8 subjects for BOLD activated region) is consistent with the increased BOLD contrast at 7 T and a smaller voxel volume causing less partial voluming. The 'combined mask' giving the region activated in both BOLD and CBF gave a lower value (M = 28 ± 2 %, 4 subjects), which would be expected due to the exclusion of large draining veins from this region.
All previous calibrated BOLD experiments have assumed a coupling relationship between CBV and CBF in Eq. (1), most using the relationship that was originally measured in rhesus monkeys, using PET and hypercapnia (α = 0.38) (Grubb et al., 1974). However since deoxygenated haemoglobin is the source of BOLD contrast, the venous, rather than total blood volume compartment should be used in the calibrated BOLD model. The venous compartment reacts passively to flow changes, whereas the arterial compartment actively constricts and dilates and so the coupling relationships for these two compartments would be expected to be different. Furthermore the global vascular response to a hypercapnic challenge may be different to the local response to neuronal activity, and the coupling relationship could change in pathology. Despite recent improvements in the understanding of this coupling (Lee et al., 2001;Kim et al., 2007;Pike, 2009, 2010), it remains preferable to avoid use of an assumed coupling constant.
Estimating q h from P ET O 2 requires an assumed value of Q 0 which is equivalent to OEF (0.4 used here; see Eq. (10) and central columns on Table 3). OEF is thought to be fairly constant across healthy subjects, but it may vary in patients. The estimate of OEF does not affect the estimated change in absolute venous blood oxygenation, ΔQ h , but it will bias the relative change q h , since this is normalised to OEF. The calculation of rvCBV is independent of errors in OEF, but a systematic error in q act will occur, which will be linearly proportional to the error in OEF, and this problem is common to all hyperoxia based methods of estimating rCMRO 2 that use a literature value for OEF. An alternative approach is to estimate global OEF from venous blood oxygenation (Lu and Ge, 2008;Jain et al., 2010). We attempted this post hoc analysis here, using the change in blood susceptibility on hyperoxia to provide a method of estimating the baseline susceptibility in the vessel independent of any model of the pattern of the field shift around the vessel. The variability in measured Q 0 (equivalent to OEF) across subjects from this technique is most probably dominated by errors in the technique, such as line profile selection and the underlying quality of the phase image, rather than biological inter-subject variability. Higher resolution gradient echo images and alternative data analysis approaches such as 2D rather than 1D profiles around the vessel, are likely to provide more reliable estimates of Q 0 . Although OEF is thought to be homogeneous across a healthy brain, a global OEF measurement may not be appropriate in some vascular pathologies (such as stroke), where local differences in OEF may be present, but with improved quality phase maps it will be possible to measure OEF more locally to the site of activation.
This work used a linear dependence of R 2 * on venous deoxygenation (β =1) (Yablonskiy and Haacke, 1994) which assumes that the system is in the static dephasing regime (Kennan et al., 1994). This assumption is considered reasonable at 7 T, where larger frequency shifts occur around venules compared to at lower field strengths (Boxerman et al., 1995). The linear correlation observed between R 2 * and P ET CO 2 and R 2 * and P ET O 2 at 7 T (Driver et al., 2010;Croal et al., 2012a) is consistent with this assumption. This is an extravascular signal model, which is reasonable since the short T 2 * of venous blood at 7 T (Blockley et al., 2008) means that the intravascular GE signal contribution is small at TE = 25 ms (Duong et al., 2003). At the clinically accessible field strengths of 1.5 T and 3 T, the intravascular contribution to BOLD signal becomes significant and the extravascular BOLD signal includes a significant diffusive component, so the relationship between R 2 * and Q is no longer linear (Ogawa et al., 1993;Kennan et al., 1994;Boxerman et al., 1995). An extra compartment could be added to the BOLD signal model to account for the intravascular signal component, or this could be suppressed by using bipolar diffusion gradients. The supra-linear relationship, relating R 2 * and Q by the power β, widely used for low field calibrated BOLD (Davis et al., 1998) could be adopted here (e.g. β = 1.5 at 1.5 T).
Common to any hyperoxia-based calibration study, a possible source of error is a change in arterial blood oxygenation on hyperoxia, which will contribute to the hyperoxia BOLD signal. Considering the changes in arterial and venous saturations on hyperoxia illustrated in Fig. 1b, the effect of hyperoxia on venous blood susceptibility is about five times greater than the effect on arterial blood susceptibility. Specifically an increase in arterial oxygen saturation of 0.017 (for oxygen partial pressure PaO 2 changing from 110 to 500 mm Hg) will cause a decrease in volume susceptibility (less paramagnetic, smaller shift relative to tissue) of Δχ v =−0.0018×10 −6 (cgs units), based on Δχ v (fully deoxygenated− fully oxygenated haemoglobin)=0.264×10 −6 (Spees et al., 2001) and a haematocrit of Hct=0.4. When this is combined with the fact that the venous blood volume is about three times greater than the arterial blood volume (Lee et al., 2001) this leads to approximately a 15 fold larger effect of hyperoxia on venous blood signal than on arterial blood signal. Assuming there is no change in arterial saturation on activation, then this will lead to a small constant overestimation of both M and M′ which will lead to a small underestimation of rvCBV and hence q act .
In contrast, hyperoxia will also lead to an increase in dissolved oxygen in arterial plasma which will cause an increase in arterial susceptibility (more paramagnetic, bigger shift relative to tissue). With the amount of dissolved oxygen determined by ε ⋅ PaO 2 in Eq. (9) the volume susceptibility contribution of oxygen dissolved in blood plasma is given by: where ε =0.0031 ml/(dl blood ·mm Hg) is the solubility coefficient of oxygen in blood, χ m =+3415× 10 −6 cm 3 mol −1 (CRC, 2010) is the molar susceptibility of oxygen (cgs units) and V m = 24.5 L/mol is the molar volume, the volume occupied by 1 mol of ideal gas at room temperature and atmospheric pressure. Therefore the contribution of  dissolved oxygen to blood susceptibility is χ v (PaO 2 = 110 mm Hg) = 0.0005× 10 −6 at normoxia and χ v (PaO 2 =500 mm Hg) =0.0022 × 10 −6 at hyperoxia (cgs units). The resulting increase in susceptibility of Δχ v = +0.0017 ×10 −6 will cause a small overestimation of rvCBV and q act , partially cancelling the underestimation of rvCBV and q act predicted in the previous paragraph due to increased arterial haemoglobin saturation with hyperoxia. To put both these changes in arterial susceptibility in context, the estimated change in venous oxygenation saturation of 0.068± 0.003 due to hyperoxia and 0.14 ±0.01 due to the motor task will correspond to a susceptibility change of 0.0072× 10 −6 and 0.0148× 10 −6 respectively (cgs units).
If the arterial blood volume changes on activation, this will also cause a change in signal on activation that is unaffected by hyperoxia, leading to an error in the estimate of q act (but not in rvCBV). Considering the relative changes in the arterial and venous blood volumes, and the T 2 * of arterial and venous blood and tissue, it is estimated that this will have an effect of about 5% on the difference in the intercepts in Fig. 1a that is used to estimate q act .
The model used here assumes that q act is the same at both normoxia and hyperoxia which in turn assumes negligible changes in both CBF and CMRO 2 on hyperoxia. The effect of hyperoxia on vasoconstriction and CBF is a matter of debate in the literature. Studies measuring CBF during hyperoxia have used fixed inspired gas mixtures to induce hyperoxia, and showed a decrease in CBF with hyperoxia (Kety and Schmidt, 1948;Watson et al., 2000;Kolbitsch et al., 2002;Bulte et al., 2007b). Attempts have been made to correct for this CBF decrease based on a look up table . However, as well as hyperoxia, these fixed inspired gas mixtures cause hypocapnia (reduced P ET CO 2 ), which will cause a decrease in CBF. Graded hypercapnia has been used to try to separate hypocapnic from hyperoxic effects on CBF, measured using continuous-ASL (Floyd et al., 2003), finding a decrease in CBF with hyperoxia. However, in more recent work (Zaharchuk et al., 2008), the apparent decrease in CBF measured by continuous-ASL was mostly accounted for by a change in arterial blood T 1 due to hyperoxia, rather than an actual CBF decrease. Recent work has found no change in global CBF on hyperoxia, when maintaining isocapnia during hyperoxia and measuring flow using phase contrast MRI (which is insensitive to T 1 changes) and arterial CBV (Croal et al., 2012b). The effect of hyperoxia on CMRO 2 has not been addressed in the literature and needs further investigation.
In this experiment, the change in P ET CO 2 during hyperoxia was − 0.5± 0.1 mm Hg, less than that for an equivalent fixed inspired hyperoxia method (ΔP ET CO 2 ≈−3 mm Hg during a 60% O 2 challenge ). The effect on CBF and resulting BOLD signal is expected to be insignificant (~2% decrease for CBF (Noth et al., 2006) and~0.2% decrease for BOLD (Driver et al., 2010)). However, in this study a brief overshoot in P ET CO 2 of 1-2 mm Hg occurred during the transition to hyperoxia. In most cases this overshoot recovered back to baseline by the time of the motor task, but it did overlap with some motor trials, which were therefore discarded. It is likely that the RespirAct™ algorithm could be modified to minimise this effect. However, any method inducing hyperoxia whilst maintaining isocapnia would be suitable for this technique.
The susceptibility difference at the air-tissue interface in the oral cavity and frontal sinus will increase on hyperoxia, leading to increased field inhomogeneity across the brain (Blockley et al., 2009;Pilkinton et al., 2011), which will be worse at higher field strengths. Changing the oxygen fraction of air from 21 to 60% O 2 changes the volume susceptibility of air by Δχ v =0.054×10 −6 (cgs units, assuming χ m (O 2 )=+3415× 10 −6 cm 3 mol −1 and χ m (N 2 )=−12×10 −6 cm 3 mol −1 (CRC, 2010)). For the EPI acquisition used here, the resulting change in field inhomogeneity has been shown to cause image distortion and intra-voxel dephasing close to the frontal sinus, but to have small effects in the motor cortex that would not significantly affect the results presented here (Driver et al., 2011).

Conclusions
A new approach to BOLD calibration has been proposed and implemented, where a task is performed both at normoxia and hyperoxia. This uses hyperoxic BOLD contrast to estimate the change in venous blood oxygenation and venous CBV during a motor task. OEF was estimated from the phase data, allowing the effect of hyperoxia on venous blood oxygenation to be estimated on a subject by subject basis. By including an ASL measurement, the relative change in CMRO 2 was estimated. The measured change in CMRO 2 agrees with previous work, using other methods, but unlike those Fig. 6. Voxelwise maps of M, rvCBV and q act across the combined mask (the intersection of BOLD and CBF motor activation masks). methods, this approach does not assume a poorly characterised relationship between vCBV and CBF.