Field strength dependence of grey matter R 2 * on venous oxygenation

The relationship between venous blood oxygenation and change in transverse relaxation rate (ΔR2*) plays a key role in calibrated BOLD fMRI. This relationship, defined by the parameter β, has previously been determined using theoretical simulations and experimental measures. However, these earlier studies have been confounded by the change in venous cerebral blood volume (CBV) in response to functional tasks. This study used a double-echo gradient echo EPI scheme in conjunction with a graded isocapnic hyperoxic sequence to assess quantitatively the relationship between the fractional venous blood oxygenation (1-Yv) and transverse relaxation rate of grey matter (ΔR2⁢GM*), without inducing a change in vCBV. The results demonstrate that the relationship between ΔR2* and fractional venous oxygenation at all magnet field strengths studied was adequately described by a linear relationship. The gradient of this relationship did not increase monotonically with field strength, which may be attributed to the relative contributions of intravascular and extravascular signals which will vary with both field strength and blood oxygenation.


Introduction
The blood oxygenation level dependent (BOLD) signal arises from the interplay between cerebral blood flow (CBF) and the cerebral metabolic rate of oxygen (CMRO 2 ) consumption. However, the magnitude of the BOLD signal also depends on baseline physiological parameters such as cerebral blood volume (CBV), oxygen extraction fraction (OEF) and haematocrit (Hct) (Kim et al., 2012;Levin et al., 2001;Lin et al., 1998). These parameters can vary for a given individual across repeated scan sessions, between healthy subjects, and in particular across patient groups. This results in variability in the magnitude of the BOLD signal, which can be minimised by using the calibrated BOLD signal to measure changes in the CMRO 2 non-invasively. The original calibrated BOLD model of Davis et al. (1998) provides a non-invasive method to quantify the fractional change in cerebral metabolic rate of oxygen (CMRO 2 ) by using the vasodilatory effect of hypercapnia to normalise the BOLD signal.
Current models of calibrated BOLD depend on the relationship between the apparent transverse relaxation rate (R * 2 ) and concentration of deoxyhaemoglobin (dHb) in the blood, or fractional venous oxygenation (1−Y v ) (Ogawa et al., 1993). The contribution to R 2 * from deoxygenated blood (Δ R * 2 ) can be described as a function of venous CBV (vCBV) and blood concentration of dHb, where [dHb] is a function of OEF and Hct (Davis et al., 1998), and k is a constant incorporating properties such as vessel geometry, magnet field strength (B 0 ) and blood-tissue susceptibility differences. The parameter β has traditionally been predicted from Monte Carlo simulations and analytical models (Boxerman et al., 1995;He et al., 2007;Ogawa et al., 1993;Yablonskiy and Haacke, 1994). This involves modelling the microscopic magnetic field gradients between deoxygenated blood vessels and the surrounding diamagnetic tissue to estimate the resulting enhancement of the R * 2 decay rate. For the extravascular (EV) compartment, the result depends on whether the system is in the static dephasing regime (larger vessels, β →1) or the diffusive dephasing regime (smaller vessels, β→2). The vessel size limit where the transition between these regimes occurs reduces as B 0 increases (Kennan et al., 1994), such that Monte Carlo simulations predict the widely accepted value of β=1.5 for a vascular network at 1.5 T, and β=1 at 7 T (Boxerman et al., 1995;Ogawa et al., 1993;Yablonskiy and Haacke, 1994). However, the reliability of these results depends on how realistic the model is. In particular, models should also consider the effect of changes in intravascular (IV) signal contributions due to concomitant changes in vCBV, for example due to vessel dilation in response to functional activation or hypercapnia, and vessel constriction in response to non-isocapnic hyperoxia. This IV contribution can be significant as the R * 2 of blood increases quadratically with increasing deoxygenation (Blockley et al., 2012;Li et al., 1998). At fields below 3 T the intravascular signal from the veins represents a significant fraction of the BOLD signal, estimated to be approximately 57% at 1.5 T, 36% at 3 T, and negligible at 7 T (Uludag et al., 2009;van der Zwaag et al., 2009), depending on echo time.
More recently, experimental measurements of β have been published Rossi et al., 2012;van der Zwaag et al., 2009;Wise et al., 2013), which are in reasonable agreement with model-based predictions. However, while theoretical estimates may not fully account for intravascular (IV) signal contributions, experimental measurements to date have been subject to potential concomitant changes in vCBV, arising from vessel dilation in response to functional activation or hypercapnia, or vessel constriction in response to nonisocapnic hyperoxia.
This study aims to measure the relationship between fractional change in concentration of deoxyhaemoglobin in venous blood and grey matter R * 2 , without a confounding change in vCBV, using a dynamic isocapnic hyperoxic stimulus at 1.5, 3 and 7 T.

Material and methods
The study was approved by the University of Nottingham Medical School Ethics Committee and all subjects gave prior informed written consent. Eight healthy volunteers (aged 24-27 years, 5 male) participated in the study. MR data was collected on Philips Achieva systems at three B 0 s, a 1.5 T system with 8-channel SENSE receive coil, a 3 T system with 32-channel SENSE receive coil and a 7 T system with head volume transmit and 32-channel SENSE receive coil.
An isocapnic hyperoxic stimulus was delivered to subjects using a feed-forward, low gas flow system (RespirAct™, Thornhill Research Inc., Toronto, Canada) in which arterial blood gases were controlled by a prospective targeting system consisting of a computerized gas blender and sequential gas delivery breathing circuit. The system uses the method of Slessarev et al. (2007) to target end-tidal PCO 2 (P ET CO 2 ) and PO 2 (P ET O 2 ) independently of each other and independent of breathing pattern (Fierstra et al., 2013). Subjects received medical air via the RespirAct™ until the respiratory challenge commenced. The respiratory paradigm consisted of three minutes of normocapnic (maintaining P ET CO 2 at resting levels) normoxia, followed by a 4 minute normocapnic linear increase in P ET O 2 to a maximum level of 500 mmHg (an approximate 8 mmHg increase in P ET O 2 with each breath), a linear decrease in P ET O 2 to return to baseline over 4 minutes while continuing to maintain normocapnia, and a final 1 minute of normocapnic normoxia. The total respiratory paradigm was thus 12 minutes in duration.
Double echo gradient echo EPI (GE-EPI) data were acquired throughout the respiratory paradigm. At all three magnet field strengths, five contiguous axial slices positioned superior to the ventricles were acquired, with an in-plane resolution of 2×2 mm 3 , 3 mm slice thickness, and field of view of 212×188×15 mm 3 , 212×184×15 mm 3 and 192×192×15 mm 3 at 1.5, 3 and 7 T respectively. Imaging parameters were TR of 3 s, SENSE factor 2.5, flip angle 90°, and TE 1 /TE 2 of 25/80 ms, 22/60 ms and 16/47 ms at 1.5, 3 and 7 T respectively. These double echo GE-EPI echo times were selected to optimally fit T 2 * based on the expected grey matter T 2 * value at each B 0 , as reported in Peters et al. (2007). In addition, at normoxia, a series of inversion recovery GE-EPI images were acquired at 10 inversion times (100, 200, 300, 500, 700, 1000, 1300, 1500, 2100 and 2600 ms) with matched geometry to the double echo GE-EPI data, which were fit for longitudinal relaxation time, T 1 . The first echo (TE 1 ) images of the double echo GE-EPI datasets were motion corrected using MCFLIRT (FSL, fMRIB, Oxford, UK) and the resulting transforms were applied to the second echo (TE 2 ) images. A grey matter (GM) mask was formed by thresholding T 1 maps at 1.0≤T 1 ≤1.3 s, 1.45≤T 1 ≤1.75 s, and 1.8≤T 1 ≤2.1 s for 1.5, 3 and 7 T respectively. To remove any remaining contributions in the GM mask from cerebral spinal fluid (CSF) and large veins, those voxels with highest and lowest 3% of signal intensity respectively in the second echo (TE 2 ) GE-EPI images were removed from the GM mask. R * 2 maps were then calculated within the GM mask on a voxel-wise basis using where S 1 and S 2 denote the signals measured at TE 1 and TE 2 . The mean baseline R * 2 was then estimated for the GM mask from the normoxic period, and the baseline R * 2 subtracted from each R * 2 value on hyperoxia, to obtain a timeseries of ΔR * GM 2 . The P ET O 2 trace was linearly interpolated to a TR of 3 s and P ET O 2 measurements were converted to venous oxygenation (Y v ), assuming a baseline OEF of 0.4 (Chiarelli et al., 2007). Y v values were binned across the full theoretical range of 0.6 to 0.7, in steps of 0.025 to enable the corresponding R * GM 2 data to be group averaged over these bins. For each individual subjects' data, at each B 0 , the average value of ΔR * GM 2 within the GM mask was fitted to two models. First, a simple linear function (1), and A and C are the fitted gradient and intercept (which takes account of the fact that we do not measure ΔR 2 * GM with respect to Y v =0). Second, to a power law function to estimate A′, β and C′, using a non-linear, three parameter least squares fit. Data were also fitted to Eq. (4), fixing β to published B 0corrected values (1.5, 1.3 and 1 for 1.5 T, 3 T and 7 T respectively (Bulte et al., 2009;Davis et al., 1998;Driver et al., 2012)). The Chisquared (Χ 2 ) was used as a measure of the goodness of fit.
We also simulated the relationship between the EV signal and venous oxygenation by using the known blood deoxyhaemoglobin relaxivity to separate the IV component from the total tissue signal, to explore contributions of any residual IV signal which was not accounted for by removal of voxels containing large vessels. While the intravascular relaxation rate R * IV  Table 1 Average baseline and peak P ET O 2 values and corresponding venous oxygenation (Y v ) estimated according to the interpolation method outlined in Chiarelli et al. (2007). Baseline R S * ( ) GM 2 −1 and the gradients (A) of ΔR * GM 2 versus (1−Y v ) for both β=1 (linear model) and typical B 0 -adjusted values of β. Χ 2 depicts goodness-of-fit. Data (mean ± SEM) averaged across subjects for each magnet field strength (B 0 ).

B 0 (T) P ET O 2 (mmHg)
Venous oxygenation (Y v ) MR parameters for GM mask Χ 2 for nonlinear fit using B 0 -corrected β values (Eq. (4)) Χ 2 for nonlinear 3 parameter fit (Eq. (4) Assuming that ΔR * EV 2 also varied linearly with Y over the range of interest and a literature value of f, the variation in S 1 /S 2 versus (1−Y v ) . In this simulation the differences in T 1 and spin density between tissue and blood were assumed negligible; in fact tissue possesses a slightly lower spin density and shorter T 1 than blood so these effects will work in opposing directions.

Results
Average baseline and peak P ET O 2 values achieved at each B 0 are shown in Table 1, along with the corresponding Y v value. Group averaged P ET O 2 traces are shown in Fig. 1A-C. P ET CO 2 values were maintained within ± 0.39 mmHg of an individuals' resting level, with the mean change between baseline normoxia and peak hyperoxia being 0.70 ± 0.12 mmHg, across all subjects and B 0 . Such fluctuations in P ET CO 2 were assumed to have minimal influence on R * 2 . Baseline (normoxic) R * GM 2 increased significantly with B 0 (p < 0.005, Wilcoxon signed rank test). Example Y v and ΔR * GM 2 timecourses are shown for representative subjects across B 0 in Fig. 1D Table 1.  NeuroImage 146 (2017) 327-332 This was significantly larger for 7 T than 1.5 or 3 T (p=0.043, Kruskall-Wallis one-way ANOVA), although no significant difference was observed between 1.5 and 3 T. The data was sufficiently described by the linear 2 parameter (A and C) model fit (Eq. (3)) (p < 0.01 in all cases). The 3 parameter model, showed no significant improvements in goodness-of-fit (Χ 2 Difference test) when compared to either of the 2 parameter fits (linear, and fixed B 0 -corrected values of β).
When R * EV 2 was estimated across binned Y v values based on the measured signals and estimated R * IV 2 , the results demonstrated a monotonic increase in the gradient of ΔR * EV 2 and (1−Y v ) with B 0 (Fig. 3B). Fig. 4 illustrates subsequent intravascular and extravascular GM signal components for the echo times and field strengths used in this study, for both normoxia and hyperoxia, using Eq. (5).

Discussion
The relationship between ΔR * GM 2 and fractional venous oxygenation (1−Y v ) was measured at 1.5, 3 and 7 T, and found to be welldescribed by a linear relationship at all B 0 s. This contrasts with previous studies which inferred that β≥1 by manipulating B 0 Rossi et al., 2012;van der Zwaag et al., 2009), and by the direct measurement of a value of β of 1.35 at 3 T using the response to combined hyperoxia and hypercapnia (Wise et al., 2013). In our study β was measured using an isocapnic hyperoxic stimulus and so the results are not confounded by changes in vCBV which will have occurred in previous studies, with vCBV altered by motor activation (van der Zwaag et al., 2009), the vasodilative effect of hypercapnia Wise et al., 2013) or non-isocapnic hyperoxia (Rossi et al., 2012).
Baseline R * GM 2 increased monotonically and approximately linearly with B 0 , in agreement with published cross-field values (Rossi et al., 2012;van der Zwaag et al., 2009). In contrast, the gradient of ΔR * GM 2 versus (1−Y v ) measured in this study did not increase monotonically (Fig. 3). Extravascular dephasing depends on the intravascular frequency shift (Δω), which depends linearly on both (1−Y v ) and B 0 : where k ″ is a constant reflecting vessel geometry, γ is the gyromagnetic ratio and Δχ is the difference in magnetic susceptibility between fully deoxygenated blood and tissue. Therefore, if we were only sensitive to extravascular signal in the static or indeed diffusive regimes, the gradient of ΔR * GM 2 versus (1−Y v ) would be expected to increase monotonically with B 0 (Boxerman et al., 1995;Ogawa et al., 1993;Yablonskiy and Haacke, 1994). In practice, the gradient was effectively constant between 1.5 and 3 T and increased faster than linearly with B 0 between 3 and 7 T (Table 1).
In practice, the non-negligible intravascular (IV) signal will contribute to the signal from GM (Hare et al., 2016), even when signal from large vessels has been masked out, as was the case in this work. Therefore, the unexpected variation in the gradient of ΔR * GM 2 versus (1−Y v ) with B 0 might be explained by a variation in the ratio of IV signal to EV signal contributions with both B 0 and echo time (TE). Here we combined measured values of ΔR * GM 2 with previously published values of R * IV 2 to estimate ΔR * EV 2 , and the results indicated that the gradient of ΔR * EV 2 versus (1−Y v ) increased monotonically with B 0 as expected (Boxerman et al., 1995;Ogawa et al., 1993;Yablonskiy and Haacke, 1994).
If a multi-compartment signal decay is fitted to a single exponential, the apparent R * 2 will depend on the choice of echo times at which the signal is sampled. Fig. 4 illustrates how the estimated relative contribution of intravascular signal to total signal varies with oxygenation and field strength for the T 2 * values measured here. The IV signal is smaller at normoxia than at hyperoxia, and therefore contributes less to the total signal as (1−Y v ) increases. Furthermore, for the echo times used here, the IV signal is more sensitive to changes in fractional venous oxygenation than the EV signal, particularly at lower field. This varying intravascaular signal contribution may explain the unexpected variation of the gradient of ΔR * GM 2 versus (1−Y v ) with B 0 . However, it should also be noted that any trends in the gradient of a linear fit over a limited range of (1−Y v ) could be hard to interpret if the underlying function was not in fact linear. Future studies designed specifically to explore these relationships should use a multi-echo EPI readout and vascular crushing to separate the different IV and EV compartments. In this study the frequency offset (Δω) between tissue and blood was varied by changing both venous oxygenation and B 0 . To fully characterise the dependency of grey matter R * 2 on Δω, a larger intravascular susceptibility perturbation is required, either by modulating venous blood oxygenation over a wider range, which is difficult to achieve physiologically, or by using a paramagnetic contrast agent.
In summary if the underlying physical value of β for any compartment is non-linear, then the effective value of β will depend on experimental parameters and local tissue microvasculature (Griffeth et al., 2011;Martindale et al., 2008), but for calibrated BOLD an additional experiment to determine β would be impractical. However for the experimental protocol used here, which was deliberately characteristic of standard fMRI experiments, a linear approximation was sufficient to describe the observed behaviour across the limited range of (1−Y v )=0.3-0.4 relevant to the fMRI BOLD effect at all field strengths, which could simplify the calibration of the BOLD signal (Blockley et al., 2015;Griffeth et al., 2013;Merola et al., 2016). Recent work (Griffeth et al., 2011;Merola et al., 2016) has tested calibrated BOLD models by comparing the relatively simple BOLD signal models with a more detailed, multi-compartment BOLD signal model, treating this more detailed model as a gold standard from which to test the specificity of the calibrated BOLD approach to simulated data. The main conclusion of this work was that the physical meaning of the parameter β (alongside a physiological CBV/CBF coupling parameter) should be relaxed and their values treated as free parameters.
It has previously been shown that a strong hyperoxic stimulus (FiO 2 =1) will cause macroscopic alterations in inferior brain regions at 3 T which can increase the measured R * 2 , although such effects were reduced under moderate hyperoxia (FiO 2 =0.5) (Pilkinton et al., 2011). Similarly at 7 T, we have previously observed whole-head field shifts under moderate hyperoxia (~FiO 2 =0.6, P ET O 2 =500 mmHg), although this only resulted in significant changes in R * 2 in regions close to the frontal sinus (Driver et al., 2011). In this work such macroscopic effects were mitigated by selecting slices superior to the ventricles, and using moderate hyperoxia. Nonetheless, further quantification of macroscopic inhomogeneities caused by hyperoxia and their impact on R * 2 are required to enable whole-head analysis. Furthermore, a different paradigm (e.g. sinusoidal ) would make it easier to remove any drifts from the data.

Conclusions
This work presents the first experimental measurement of β in humans using an isocapnic hyperoxic stimulus to remove confounds associated with a change in CBV arising from functional activation or modulation of P ET O 2 during a respiratory challenge. It was also found that for the regimes relevant to calibrated BOLD experiments, a value of β=1 is a reasonable assumption, which may help simplify the calibrated BOLD model.