A General Model to Calculate the Spin–Lattice Relaxation Rate (R1) of Blood, Accounting for Hematocrit, Oxygen Saturation, Oxygen Partial Pressure, and Magnetic Field Strength Under Hyperoxic Conditions

Under normal physiological conditions, the spin–lattice relaxation rate (R1) in blood is influenced by many factors, including hematocrit, field strength, and the paramagnetic effects of deoxyhemoglobin and dissolved oxygen. In addition, techniques such as oxygen‐enhanced magnetic resonance imaging (MRI) require high fractions of inspired oxygen to induce hyperoxia, which complicates the R1 signal further. A quantitative model relating total blood oxygen content to R1 could help explain these effects.

of R1 in blood is important for several research techniques such as arterial spin labeling and black blood imaging. 1,2 Under normal physiological conditions, the relaxation rate in blood is influenced by many factors, including hematocrit, field strength, temperature, protein concentrations, and the weak paramagnetic effects of deoxyhemoglobin and dissolved oxygen in the plasma. [3][4][5] In addition, evolving techniques such as oxygen-enhanced magnetic resonance imaging (MRI) require high fractions of inspired oxygen to induce hyperoxia, which complicates the R1 signal further. [6][7][8][9][10][11] Whilst a linear increase in R1 with oxygen partial pressure has been reported in phantom experiments, measurements in blood have produced contradictory effects, in some cases showing an unexpected decrease in R1 following induction of hyperoxia. 12 In response to these contradictory results, a recent Journal of Magnetic Resonance Imaging Editorial by Desmond et al suggests that "future work could expand the interpretation of the T1 and T2* changes using quantitative models relating total blood oxygen content, partial pressure of oxygen (PO 2 ) and oxygen saturation (SO 2 ) to the MRI measures." 13 Existing models for predicting R1 in blood can account for SO 2 , hematocrit, and magnetic field strength, but are not applicable under conditions of artificially increased PO 2 , or hyperoxic experiments. 14 Therefore, in this paper our aim was to establish a general model that could estimate the R1 of blood, accounting for hematocrit, SO 2 , PO 2 , and magnetic field strength under both normal physiological and hyperoxic conditions.

Background Theory
The relationship between longitudinal relaxation and the paramagnetism of oxygen in blood has a long history in MRI, beginning when Young et al 15 and Tripathi et al 16 found that increasing inspired oxygen fractions led to a small but definite increase in the R1 of blood in the left ventricle (i.e., oxygenated blood), but no change in R1 in the right ventricle (i.e., deoxygenated blood). To identify the physiology underlying these contrast changes, experiments using the vitreous fluid of the eye 17 and ex vivo blood samples at varying levels of oxygenation led to the hypothesis that the paramagnetic dissolved oxygen in the plasma provides the contrast during hyperoxia. 18 This hypothesis was confirmed in ex vivo blood samples by d'Othée et al, 5 where measurements of the bound and dissolved oxygen plotted alongside the measured R1 (reproduced in Fig. S1 in the Supplemental Material) showed that following full oxygen saturation, the R1 of blood increases linearly with increasing dissolved oxygen, while the bound oxygen remains constant. Prior to full oxygen saturation, however, experiments by Silvennoinen et al showed that the main source of contrast change is caused by changing levels of paramagnetic deoxyhemoglobin, where R1 decreased linearly as SO 2 increased (i.e., with decreasing deoxyhemoglobin concentration). 4 In 2016, Hales et al published a general model to calculate the spin-lattice relaxation time constant (T1) of blood, accounting for hematocrit, SO 2 , and magnetic field strength. 14 The model built upon a two-compartment model proposed by Grgac et al, in which a fast-exchange model 19 model was used to describe blood R1 because the water exchange rate between the plasma and cytoplasm is faster than the difference in R1 between these two compartments, [20][21][22] given by: where R1 b is the relaxation rate of whole blood, SO 2 is the oxygen saturation fraction (range of 0-1), R1 p is the longitudinal relaxation rate of plasma (s À1 ), and R1 e is the longitudinal relaxation rate of erythrocytes (s À1 ).
The variable f e is the fraction of water in whole blood that resides in erythrocytes (0-1), which is described by the equation: where Hct is the hematocrit (range of 0-1), W RBC is the volume fraction of water within the erythrocyte (typically valued at 0.70 due to hemoglobin occupying approximately 30% of the erythrocyte volume), and W plasma is the volume fraction of water within the plasma (typically valued at 0.95, where the other 5% volume is occupied by plasma proteins such as albumin). 21 Hales et al extended this model to account for the effect that deoxyhemoglobin has on R1 e , as a paramagnetic contrast agent, expressing R1 e as a function of SO 2 : where R1 eox is the relaxation rate of erythrocytes when SO 2 = 100%, [Hb] is the mean corpuscular hemoglobin concentration (5.15 mmol Hb tetramer/L plasma), and r1 dHb is the molar relaxivity of deoxyhemoglobin (in s À1 L plasma in erythrocyte/mmol Hb tetramer). Combining Eqs. 1-3 results in the following model for calculating R1 of blood (Eq. 4): In addition, Hales et al modeled R1 eox , R1 p , and r1 dHb to have a linear dependence on B0 (where β 0 and β 1 are the y-intercept and slope for the linear fit): To address this limitation, we defined a new parameter, r1 pOx , the relaxivity of dissolved oxygen in the plasma in s À1 /mmHg oxygen (Eq. 8), which we modeled as having a linear dependence on field strength (Eq. 9) in the same manner as Hales et al had for r1 dHb , R1 eox , and R1 p , (Eqs. 5-7): Although we model Eqs. 5-7 and 9 as linear, the true relationship between field strength and both R1 and the relaxivity (r1) and is nonlinear and depends on the material. [23][24][25] For R1, protons in highly mobile molecules such as water or plasma will be less affected by changing field strength, whereas in materials such as biological tissues the R1 has been empirically shown to decrease as field strength increases by approximately B0 À 1 /3 (see Fig. S2a in the Supplemental Material). 23 However, when examined from the range of 1.5-7 T, measurements by Rooney et al (provided in Fig. S2b in the Supplemental Material) suggest that using a linear model for R1 eox (B0) (Eq. 5) and R1 p (B0) (Eq. 6) is appropriate in the field strength range of this model (1.5-8.45 T). Similarly, an experiment by Chou et al measured the relaxivity (r1) of paramagnetic contrast agents over a large range of field strengths (0-12 T), revealing a Lorentzian shape (reproduced in Fig. S3a in the Supplemental Material). 24 However, when only examined between the range of 1.5-8.45 T, the measurements by Chou et al (provided in Fig. S3b in the Supplemental Material) suggest that using a linear model for r1 dHb (B0) (Eq. 7) and r1 pOx (B0) (Eq. 9) is appropriate in the field strength range used in this model (1.5-8.45 T). Therefore, a linear model has been chosen for Eq. 5-7 and 9. Equation 8 is a function of both SO 2 and PO 2 , which can be related by the invertible form of the hemoglobin-oxygen dissociation curve (the Hill equation) 26 : where P50 is the half-saturation PO 2 and n is the Hill exponent for hemoglobin (typically 2.7). 27 P50 will vary by species and can be set to match the species of interest by using values obtained experimentally. 26 By substituting Eq. 10 into Eq. 8, we obtain Eq. 11 for calculating the R1 of blood that is no longer limited to normal physiological oxygen ranges. Thus, Eq. 11 can be used to calculate R1 even when the partial pressure of oxygen in the plasma reaches normobaric hyperoxic levels.
In addition, it is known that many paramagnetic contrast agents interact with macromolecules and exhibit a higher relaxivity at higher macromolecule concentrations. 28 This is also true for molecular oxygen, and blood with artificially high levels of hematocrit (Hct = 0.6) has been shown to have a higher relaxivity than blood at regular hematocrit levels. 29 To account for this hematocrit sensitivity, we also tested the effect of adding an extra term, β 1,r1pOxHct , to Eq. 9, assuming a linear dependence on hematocrit.

Data Analysis
In order to estimate values for the parameters in the new model, 126 data points for R1 in blood were collected from the literature with measured or estimated Hct, SO 2, and PO 2 at a range of field strengths (1.5-9 T) and under normoxic and hyperoxic conditions (Table 1). 3,5,20,[29][30][31][32][33][34] Data were extracted either from tables, manuscript text or graphs, and if data extraction from graphs was necessary, a digital plot analyzer was used to reliably extract values. 35 The distribution of the resulting data is shown in Fig. S4 in the Supplemental Material. An accurate measurement of T1 depends on the selection of acquisition protocols and parameters (such as repetition and inversion times), and therefore the acquisition details from each experiment are listed in Table S1 in the Supplemental Material.
The T1 values from Stefanovic and Pike and d'Othée et al were increased by 12% to account for the shift in T1 between 22 C and 37 C. 10,35 If not reported, PO 2 values were calculated from the reported SO 2 values using Eq. 9. Since the value of P50 (halfsaturation PO 2 ) can differ greatly between species, a P50 of 25 mmHg was used for all human data points, and a P50 of 44 mmHg was used for the bovine and rat, based on the P50 values found experimentally in human, bovine, and rat blood samples. 26 For model fitting, the SciPy optimize function for nonlinear least-squares fitting was used. 36 Equation 11 was fitted using a randomized subset of 95% of the 126 data points, fitting with the R1, B0, P50, Hct, and PO 2 values simultaneously, and testing on the The distribution of the data points is plotted in Fig. S1 in the Supplemental Material. The MRI acquisition details from each experiment are listed in Table S1 in the Supplemental Material.
remaining 5%. The dataset was split for fitting and testing using the SciKitLearn train_test_fit function with shuffling. 37 This process was iterated 10,000 times, and the mean and 95% confidence interval of the distribution of fitted values for each of the parameters were recorded (listed in Table 2).

Effect of W RBC and W plasma Constants
In Eq. 2, the values used for W RBC and W plasma were 0.70 and 0.95, respectively. Experimentally, W RBC has been measured to be at maximum 78% and at minimum 54%, 38 and 94% has sometimes been used for W plasma instead of 95%. 4 To examine the effect of this range of possible variables on the final R1, the resulting R1 was calculated for all combinations of W RBC = 54-78% and W plasma = 94-96%, and Hct = 0.38-48.

Statistical Analysis
Since Eq. 11 contains eight fitting parameters, the fitting process was repeated for all combinations of fewer parameters (i.e., removing the dependence of R1 p on B0), and the Akaike Information Criterion (AIC) was calculated for each version of the model-the best-fit model according to the AIC is the model that explains the greatest amount of variation using the fewest possible independent variables. The model that scored the highest according to the AIC was therefore selected as the best model and used in further investigation. The R 2 coefficient of determination and mean squared error (MSE)

Akaike Information Criterion
The AIC score, R 2 , and MSE results of the four highest AIC-scoring models tested are listed in Table S2 in the Supplemental Material, and results are plotted in Fig. S5 in the Supplemental Material. The version of the model with β 1,R1p set to 0 (i.e., with Eq. 6 set to R1 p = β 0,R1p , hence 7 total parameters) resulted in the highest AIC score, and was therefore used in this manuscript.
The version of the model with the extra fitting parameter β 1,r1pOxHct in Eq. 9 resulted in a similar fit (R 2 = 0.927 vs. R 2 = 0.926), but the best-fit model according to the AIC was the model without this extra parameter (see Table S2 in the Supplemental Material). Therefore, the best model was the 7-parameter model which did not include this dependency on hematocrit.

Model Fitting
The modeled vs. measured R1 b values from the randomized unseen test set of each iteration are plotted in Fig. S6 in the Supplemental Material alongside the line of equality and the linear regression (R 2 = 0.91). The difference between the modeled and measured R1 b values from the randomized unseen test set of each iteration is illustrated as Bland-Altman plots in Fig. S6b-f in the Supplemental Material with the resulting mean error of R1 blood as 1 Â 10 À4 s À1 and the limits of agreement of À0.078 s À1 and 0.078 s À1 . Violin plots showing the distribution of parameter estimates from each iteration are shown in Fig. S7 in the Supplemental Material.
The median parameter values for the final optimal 7 parameters β 0,R1eox , β 1,R1eox , β 0,r1dHb , β 1,r1dHb , β 0,R1p , β 0,r1pOx , and β 1,r1pOx (with lower and upper 95% confidence intervals) are listed in Table 2. The accuracy of the final model's predicted R1 compared against the "true" R1 b values is shown in Fig. 1a. The predicted vs. true R1 b values are plotted against the line of equality (solid black line) and a linear regression (dotted line, R 2 = 0.93) in Fig. 1a, with a final MSE of 0.0013 s À2 . To examine bias in the model with respect to SO 2 , PO 2, B0, Hct, or R1 variables, Bland-Altman plots are provided in Fig. 1b-f, with the resulting mean error of R1 blood as 6 Â 10 À6 s À1 and the limits of agreement of À0.071 s À1 and 0.071 s À1 . A systematic trend in the error distribution by R1 (Fig. 1b) indicated that the model did not entirely explain the non-noise variance in the data, and possible causes of this were further examined by viewing the datapoints according to the original source (provided in Fig. S8 in the Supplemental Material).

Model Behavior
To understand the behavior of the resulting model, the effect of varying PO 2 , SO 2 , field strength, and hematocrit on R1 is illustrated using synthetic data from 0 to the upper level of normal oxygen levels (Fig. 2) and extending to hyperoxic conditions (Fig. 3).
As shown in Fig. 4, the ability of the new model to predict blood R1 was compared with the model and parameter values estimated by Hales et al, on a subset of data points (N = 128) from non-hyperoxic blood, i.e., within the range that the Hales model is stated to be valid (SO 2 = 0.4-0.99). Within this range, the model agreed well with the Hales model estimates of R1 (R 2 = 0.96), and showed slightly   The behavior of the erythrocyte and plasma components of the model is shown in Fig. 5, where the model showed that as oxygen levels increased in the blood, the R1 decreased due to the paramagnetic effect of the deoxyhemoglobin fraction, and once full saturation was reached, the predominant remaining effect was the increase in R1 due to the paramagnetic dissolved oxygen in the plasma.
The predicted changes in R1 that would be induced by an increase in oxygen levels are illustrated in Fig. 6a,b, which  shows the so-called "contradictory" negative ΔR1 that has been observed by oxygen-enhanced MRI researchers in hypoxic tumors. [6][7][8][9][10][11] As can be seen in Fig. 6c, when starting from a low saturation, the R1 first decreases until full saturation is reached, and then increases linearly with PO 2 . In Fig. 6d, the model showed that a change in PO 2 from 39 mmHg to 48 mmHg (such as in venous blood) induced a negative ΔR1, while a change (such as in aterial blood) from 98 mmHg to 600 mmHg induced a large positive ΔR1.

Effect of W RBC and W plasma Constants
The distribution of resulting R1 values for all combinations of W RBC = 54-78% and W plasma = 94-96%, and Hct = 0.38-48 is plotted in Fig. S10 in the Supplemental Material, and the resulting SD of the R1 values for a set hematocrit (range 0.38-0.48) was 0.02 s À1 , and a maximum difference of maximum difference in R1 of 0.07 s À1 .

Discussion
We have presented a general model to estimate the R1 of blood, accounting for hematocrit, oxygen saturation (SO 2 ), the partial pressure of oxygen (PO 2 ), and magnetic field strength under both normal physiological and hyperoxic conditions. The model was valid for the following parameter ranges: B0 = 1.5-8.45 T, SO 2 = 0.40-1, PO 2 = 30-700 mmHg, and adjustable for different species by adjusting the P50 constant appropriately. The model estimates agreed well with literature values (R 2 = 0.93, MSE = 0.0013 s À2 ) and the model performance was consistent across R1, SO 2 , PO 2 , B0, and Hct levels.
The new model builds upon a previous model for R1 or T1 of blood, 14 accounting for hematocrit, field strength, and SO 2 , valid between 1.5-7 T and SO 2 of 40-99% (nonhyperoxic conditions). The new model showed improved accuracy and preference within the valid range of the previous model, likely due to the four parameters added which better reflected the true behavior, and also due to a 3Â increase in data points available for model fitting. The new model was also valid over a slightly larger range of field strengths (1.5-8.45 T). The main improvement of the new model, however, was that it could account for changes in PO 2 under hyperoxic conditions.
Our model showed that there are two competing effects on blood R1 that arise from increasing oxygen levels (paramagnetic oxygen and paramagnetic deoxyhemoglobin), and that the effect on R1 due to deoxyhemoglobin dominates during lower oxygen levels. O'Connor et al reported that when breathing 100% oxygen, there is a small increase in the venous PO 2 from 39 mmHg to 48 mmHg, while the arterial PO 2 shows a high increase from 98 mmHg to 600 mmHg. 39 By modeling these hyperoxic conditions, we were able to quantify the expected ΔR1 in venous and arterial blood after breathing 100% oxygen.
By modeling venous and arterial conditions, this model also provided an explanation for the unexpected additional increase in T1 (decrease in R1) from the oxygen-enhanced drink vs. the control drink reported by Vatnehol et al in the hepatic portal vein, where the authors expected an increase in R1 from increased oxygenation. 12 The authors hypothesized that this contradictory result was due to oxygen increasing the amount of water absorbed into the bloodstream from the beverage and hence increasing the dilution effect. 12 However, according to our model, if the drink increased the oxygen saturation levels in the vein, then a decrease in R1 should be expected-which is consistent with their observations.
As mentioned in the background theory, the R1 of protons in highly mobile molecules such as plasma will be less affected by changing field strength, whereas in low mobility materials (such as protons in biological tissues) the R1 has been empirically shown to decrease as field strength increases. 23 The resulting slope of the fieldstrength dependence for R1 p (plasma) and R1 eox (erythrocytes) agreed well with this theory-the resulting slope for R1 P (i.e., β 1,R1p ) was almost 0. In fact, the version of the model with β 1,R1p set to 0 (i.e., with Eq. 6 set to R1 p = β 0,R1p , hence 7 total parameters) resulted in the highest AIC score and was therefore used in this manuscript. In addition, the removal of the dependence of R1 eox produced an extremely low AIC score, therefore suggesting that the R1 of the protons within the erythrocytes was more affected by changing field strength than the R1 of the protons within the plasma, which was consistent with this background theory.
For the purpose of this manuscript, the values used for W RBC and W plasma were 0.70 and 0.95, even though in reality these can range from W RBC = 54-78% and W plasma = 94-96%. When testing the effect of the full range of these variables on R1, it resulted in a maximum difference of 0.07 s À1 and SD of 0.02 s À1 between any combination of W RBC and W plasma values and set hematocrit, suggesting that deviation from the chosen population average values of 0.70 and 0.95 could account for the remaining variance in the R1 estimate.

Limitations
The R1 values collected were from experiments performed over a timespan of two decades with variation in the experimental equipment and T1-mapping technique used. An accurate measurement of T1 in blood samples or in vivo can be difficult and the measurement accuracy depends on the proper selection of acquisition protocols and parameters (such as repetition and inversion times). Of note, some of the R1 data points collected from Stefanovic et al were underestimated by the model-these were the only data points collected by a Look-Locker T1 mapping method, which are known to overestimate R1. 40 In addition, a systematic trend in the error distribution by R1 was noticed, however when viewing the datapoints according to the original source, it was apparent that this trend was caused due to the model consistently underestimated the values of R1 reported from one literature source, 29 in a way that was not explained by categorization by B0, Hct, SO 2 , or PO 2 values. Hueckel et al stated that the temperature of all samples was "approximately 37 C," however, inaccuracy in this measurement could have caused a higher R1 than the rest of the authors' measurements which were performed at 37 C (or corrected from room temperature).
It is important to note that the values of blood R1 have often been extracted from original plots, some hand-drawn, which is a source of error, and converted from the various original units to s À1 /mmHg, which can be another source of error. The hematocrit, SO 2 , and PO 2 measurements collected from the literature also have several limitations-eg, it was apparent that researchers tended to report certain rounded values more often than others and to different levels of precision. Altogether, these limitations are a large source of potential error in the accuracy of this model. As a future direction of this work, we have hosted the open-source model code and current R1 b measurements on GitHub (github.com/ BulteGroup/R1BloodModel) and invite the MRI community to share more modern measurements of blood R1 at a range of field strengths, hematocrit, SO 2 , and PO 2 , both in vivo and ex vivo and refit the model to improve accuracy.
The literature dataset used for the model did not include SO 2 values below 40%, so the model may not be valid below that threshold. The Bland-Altman plots showed that the prediction error was greater for a set of data points with lower SO 2 (below 50%). This could be due to the model being slightly less accurate at lower SO 2 due to not having many data points with SO 2 < 50% in the dataset (and none below 40%), or it is possible that these original reported R1 measurements were inaccurate. Regardless, we welcome future researchers to acquire new data points at lower SO 2 and refit the model to re-estimate the parameters and have more certainty within this region if it is desired.
Finally, it should be noted that there are other factors that will influence R1 in blood not included in the model. Plasma proteins that are relevant from an MRI perspective due to either diamagnetism or paramagnetism are albumin, Immunoglobulin G and fibrinogen, transferrin, and ceruloplasmin, and the concentrations of these can differ by age and sex, resulting in a varied plasma R1. 28 It is possible that some of the variability in the literature data could be accounted for if these concentrations were also measured and reported. In fact, Li et al have proposed a highly specific model of the R1 of blood accounting for hemoglobin saturation, methemoglobin, and albumin under non-hyperoxic conditions. 21 These data were not available in the literature used for the model, but if these were recorded in future studies alongside new R1 measurements they could be utilized for refitting.

Conclusions
In this work, we presented a general model to estimate the R1 of blood, accounting for hematocrit, magnetic field strength, SO 2 , and PO 2 , under both normal physiological and hyperoxic conditions, and adjustable for different species. This work was built upon a previous model, and we have both improved the accuracy of the R1 estimates within the valid range of the previous model and extended its range of validity to include hyperoxic conditions. Using the predictions of this model, we quantified and explained contradictory oxygen-induced R1 changes reported in the literature and characterized the expected R1 changes arising from the relevant physiological factors.