## Abstract

**Background** During cardiopulmonary bypass (CPB), maintaining adequate oxygen consumption (VO_{2}i) can only be achieved indirectly either by modifying oxygen delivery (DO_{2}i) through its component parts or by modulating metabolic demand through altering body temperature. The body reacts to these actions by changing OER and consequently VO_{2}i. Understanding the body’s adaptive OER dynamics can elucidate its oxygen consumption goals during CPB and help improve our ability to safely manage the patient’s journey.

**Methods** An autoregressive, integrated time-series model was trained on granular perfusion data from 879 paediatric patients (age: newborn to 18 years old) undergoing 963 CPB operations, with the outcome variable being the minute-by-minute changes in the logit transformation of OER. Variables were cardiac index, haemoglobin concentration, oxygen saturation of arterial haemoglobin and temperature. An explicit ‘disequilibrium term group’ was also included, proportional to the difference between the logarithm of VO_{2}i and logarithm of a ‘latent’ (i.e. unobserved) oxygen demand - or ‘target’ VO_{2} (tVO_{2}i) - term, with the logarithm of tVO_{2}i assumed to be a linear function of body temperature (the Van’t Hoff model). The trained time-series models were studied using permutation-based variable importance, deterministic and stochastic simulations, and subgroup analysis by acute kidney injury (AKI) grade and by temperature.

**Results** Model coefficients are consistent with an adaptive OER response to keep VO_{2}i in line with tVO_{2}i, according to body temperature. This adaptation consists of a primary rapid response for 5-10 minutes, and a secondary slow response that is estimated to last up to several hours. The model reproduces the hyperbolic shape of DO_{2}i-VO_{2}i curves - first published in 1982 - as an artefact of insufficient wait times between equilibrium-state transitions. Asymptotically, however, the model converges to a piecewise linear relationship between DO_{2}i and VO_{2}i, with supply-independence of oxygen consumption occurring above a threshold DO_{2}i. Subgroup analysis by temperature suggests that the dependence of tVO_{2}i on temperature (expressed as Q10) may be significantly stronger at low temperatures (< 28C) than at high temperatures (> 28C).

**Conclusions** This study proposes a physiologically plausible model of OER changes during CPB that is consistent with past experimental data. While during CPB, under-oxygenation is the dominant risk in the long term, slow adaptation of OER during CPB creates short-term opportunities for over-oxygenation following significant changes in variables such as cardiac index. The model provides well-defined values for tVO_{2}i at a given temperature, paving the way for further research into the effects of over- and under-oxygenation during CPB on postoperative outcomes such as AKI, and hence improvements in goal-directed perfusion protocols.

**What Is New?**

This study is the first to present a data-driven, analytical framework for predicting OER changes in response to clinical interventions during CPB.

Changes in the components of oxygen delivery cause an adaptive OER response to keep oxygen consumption in line with oxygen demand, according to body temperature.

The dependence of oxygen demand on temperature decreases as temperature increases towards normothermia, inconsistent with the accepted Van’t Hoff equation.

Children developing AKI exhibit a dampened response to changes in haemoglobin during CPB, with this dampening of response intensifying with AKI severity.

**What Are the Clinical Implications?**

This proposed, dynamic model of OER provides a novel framework for goal-directed perfusion by identifying periods of over- and under-oxygenation.

The observed, dampened response to haemoglobin changes in patients that develop AKI can be the foundation of an intraoperative tool for early diagnosis of at-risk patients.

## Introduction

During cardiopulmonary bypass (CPB), the heart-lung machine (HLM) takes over the twin functions of perfusion: oxygenation and carbon dioxide removal, and circulation of the blood. Given the central role of perfusion during surgery, much attention has been given to reducing postoperative end-organ dysfunctions via improvements in perfusion strategies, which has led to the advent of the term goal-directed perfusion (GDP). An important prerequisite for GDP is identifying an appropriate, and specific, set of ‘goals’ for perfusion. Much research has focused on the determination of a threshold for oxygen delivery - indexed to body surface area (BSA) - (DO_{2}i) for aerobic metabolism, below which the risk of acute kidney injury (AKI) increases (Ranucci *et al*, 2018; Hendrix *et al*, 2019; Dreher *et al*, 2023; Do-Nguyen *et al*, 2023; Hayward *et al*, 2023). Focusing on under-oxygenation risk is motivated by the concept of ‘supply-independence of oxygen consumption’, which posits that, as long as sufficient oxygen is delivered to the body, it is otherwise capable of extracting what it needs, without the risk of over-oxygenation. This concept is reflected in the piecewise linear (or elbow) model of oxygen consumption vs. oxygen delivery (Figure 1). We refer to the supply-independent oxygen consumption as oxygen demand, or target VO_{2}i (tVO_{2}i).

Focusing on a single threshold-DO_{2}i (tDO_{2}i) during an entire CPB session, while easy to communicate and implement by practitioners, has a few shortcomings. First, it does not take into account the temperature-dependence of oxygen demand. During CPB, especially one involving Deep Hypothermic Circulatory Arrest (DHCA), the patient’s body temperature can change by 20°C or more. Assuming a Q10 of 2.5, a 20-degree temperature change translates into a more than six-fold variation in oxygen demand during CPB. Using a single tDO_{2}i for this entire range would be either impractical (satisfying a high tDO_{2}i would clash with need to reduce flow at times for better visibility during surgery) or unsafe (allowing for a low tDO_{2}i would endanger the patient during normal temperatures), or both. Secondly, empirical evidence does not support a simple elbow model with a clear-cut tDO_{2}i. For instance, Fox *et al*, (1982) show that a hyperbolic model (Figure 1) - which is a smooth version of the elbow model - better fits the data from human and animal experiments (Paneth et al, 1957; Cheng et al, 1959; Starr, 1959)^{1}. Besides the hyperbolic function, other smooth, saturating functions such as the exponential curve have also been proposed (Lubarsky et al, 1995). Despite better fit, such smooth functions do not provide a well-defined tDO_{2}i. Ad-hoc substitutes, such as finding the intersection of the 70% confidence bands with the tVO_{2}i line (Kirklin and Blackstone, 2012), are somewhat arbitrary. One potential advantage of these experiments supporting the textbook model is that the experimental setup presumably allowed for the system to reach steady state before each measurement. However, this advantage is also the reason why such experiments are unlikely to be replicated, especially in humans. On the other hand, technological advances have allowed us to collect temporally granular perfusion data during routine CPB. In these observational settings, variables such as blood flow, haemoglobin concentration, arterial oxygen partial pressure and body temperature are all subject to constant fluctuations, driven by surgical requirements. In other words, we can never assume the system is in a steady state or equilibrium. Yet, after mathematically accounting for system dynamics, one could extract important insights from such granular data about the body and its oxygen uptake mechanisms during CPB, including the dependence of oxygen consumption on oxygen delivery, as well as the dependence of oxygen demand on body temperature. The primary aim of this study is to understand the innate adjustments to OER that occur during CPB to maintain oxygen consumption in line with oxygen demand, as dictated by the resting metabolic rate (RMR). In particular, we create a time-series model to predict the minute-by-minute changes in oxygen extraction rate (OER), as a function of the past and intended (i.e. future) changes in blood flow rate, haemoglobin level, and temperature (exogenous variables), past changes in OER (endogenous variable), and present state of the system. The secondary aim of the study was to understand how coefficients and model performance changed in patients that developed AKI compared to those that did not. Finally, we aimed to examine the relationship between oxygen demand and temperature to determine if it follows the simple parametric form of Van’t Hoff, which has been the prevailing model in the literature. This was undertaken by conducting a subgroup analysis of the data by temperature band over the entire range of temperatures recorded in our granular CPB data.

## Methods

### Study Design and Participants

This was a retrospective analysis of patients aged younger than 18 years (maximum age treated at the institution), who underwent cardiac surgery with CPB at Great Ormond Street Hospital for Children, London. The cohort was selected using a National Institute for Cardiovascular Outcomes Research (NICOR) validated data set for patients operated on between April 2019 and April 2021. All clinical data from 879 paediatric patients undergoing 963 CPB operations were collated in a research platform within the hospital’s governance structure and de-identified before analysis. Institutional approval was given for the undertaking of this project (audit number 3045). Individual consent was not required because only routinely collected de-identified hospital data were evaluated within the secure digital research environment as part of an existing research database approval (17/LO/0008).

### Data Collection

Routine clinical information was extracted from the institution’s Electronic Health Record (EHR) system using a custom structured query language script. These data included demographic information, laboratory results, intraoperative, CPB, and medication administration data, and intensive care requirements. Comorbidities (including risk scores), and outcomes were defined using Association for European Paediatric and Congenital Cardiology coding as part of the institution’s NICOR routine data submissions. Intraoperative data from the heart-lung and anaesthetic machines were captured every minute throughout surgery. Preoperative risk estimation was based on cardiac and non-cardiac preoperative risk factors using the Partial Risk Adjustment in Surgery (PRAiS2) model (Rogers et al, 2017). Baseline serum creatinine levels were taken within the 24 hours preceding surgery. Intraoperative blood biochemistry (including haemoglobin, arterial oxygen partial pressure and saturation, and mixed-venous oxygen saturation) was measured continuously using an optical fluorescence and reflectance system (CDI550, Terumo, Leuven, Belgium), calibrated every 20 minutes using a cassette-based blood gas analyzer (ABL90 Flex Plus, Radiometer, Copenhagen, Denmark). HLM variables were automatically recorded in one minute intervals via a bidirectional connection into the Institutional EHR system (Axon, Philips Capsule, Paris, France). Postoperative Acute Kidney Injury (AKI) was calculated according to the Kidney Disease Improving Global Outcomes (KDIGO) criteria (Khwaja, 2012), based on the patients’ serum creatinine change and urine levels within the first 48 hours following surgery. Patients were assigned the highest AKI grouping on the basis of these two criteria. The final KDIGO score for each patient is an integer between 0 (no AKI) and 3 (severe AKI).

### Anaesthesia

Anesthesia was induced by inhalation of sevoflurane in oxygen and, after induction, fentanyl 5 mg/kg and pancuronium100 mg/kg were given; maintenance was with isoflurane 1.0% in oxygen and air. After tracheal intubation, arterial and central venous lines were inserted. Further incremental doses of fentanyl, up to 25 mg/kg, were given during the procedure.

### CPB

The CPB circuit consisted of FX oxygenators with integrated arterial line filter and hard-shell venous reservoir (Terumo Corp). The Stöcket S5 (Stöcket, LivaNova) heart-lung machine with pole-mounted roller pumps was used with a 3T heater/chiller (Stöcket, LivaNova). The total base prime volume was 330 - 1000 mL depending on circuit size. When the predicted initial CPB hematocrit was <27% the prime consisted of packed red blood cells of the patient’s blood group (120-200 mL), a synthetic colloid (Gelofusine; B Braun Melsungen AG), heparin 1000 U/mL, 1.5 mL (Heparin Sodium; Wockhardt). Blood primed circuits were then washed with 1000 mL of balanced crystalloid solution (Plasmalyte 7.4; Baxter) by performing prebypass ultrafiltration, carried out as previously described, to ensure an initial on-CPB hematocrit of 30% (Naik and Elliott, 1993). Biochemical compatibility was then attained using sodium bicarbonate as a buffering agent. Patients requiring a clear prime received crystalloid and colloid volumes in a 40:60 mix, with 10 mL sodium bicarbonate and 2 to 5 mL heparin. All patients were systemically cooled to nasopharyngeal temperatures between 28 and 35°C unless undergoing deep hypothermic circulatory arrest, in which case the patient was cooled to 18°C. Myocardial protection was achieved with 30 mL/kg of cold blood cardioplegia (4:1 blood:cardioplegia to a final concentration of 20 mmol) of St Thomas’s Solution (IVEX Pharmaceuticals). During the rewarming phase of CPB, a gradient no greater than 10°C between the patient’s nasopharyngeal temperature probe and the heater/chiller (maximum arterial blood temperature 37.5°C) was maintained until a maximum oesophageal temperature of 36°C was achieved.

### Data Preparation

All data preparation and analysis for this research has been conducted using the R programming language and its various packages (R Core Team, 2021).

### Data Cleaning

Each operation was split into one or more ‘sessions’ of contiguous timestamps (every minute). After calculating the necessary derived fields (see below), a boolean ‘validity’ flag was calculated for each intraoperative recording, enforcing the following conditions: 0 < *Hct* < 0. 6, 0 ≤ *SvO*_{2} ≤ 100, *SvO*_{2} < *SaO*_{2}, *T >* 0, *CI >* 0, where *CI* is cardiac index (blood flow indexed to BSA) in L/min/m^{2}, *Hb* is blood haemoglobin level in g/dL, and *SaO*_{2} is the haemoglobin oxygen saturation of arterial blood, expressed as a percentage. (The last condition removes periods of DHCA, i.e., *CI* = 0, from further analysis.) Also, we flag any recordings with missingness as invalid, thus limiting our time-series analysis to valid ‘segments’ created within each session.

In calculating the change terms used for constructing the design matrix of the GARIX+ model (Equation 9), we impose a minimum segment length of 10 minutes for the GARIX+(7) model, and 25 minutes for the GARIX+(20) model. When comparing GARIX+ models of different *N*, we used a single value for the minimum segment length to prepare a common dataset.

### Derived Fields

Body surface area (BSA) is defined according to the Monsteller formula: where height is in cm, and weight is in kg.

Oxygen delivery (indexed to BSA) is given by: In the above equation, the contribution to oxygen content of blood from dissolved oxygen is ignored, thus focusing exclusively on the haemoglobin-bound oxygen. (In our data, dissolved oxygen constitutes only 1.6% of total oxygen content in arterial blood, thus supporting the plausibility of ignoring the dissolved-oxygen component in our analysis.)

Similarly, for oxygen consumption:
where, again, the contribution from change in the dissolved oxygen content of the blood is ignored. This allows us to write OER as:
Therefore, VO_{2} i can be rewritten as:

### Kirklin Equation

The following equation (page 119 of Kirklin and Blackstone, 2012) represents a hyperbolic functional form for dependence of *VO*_{2}*i* on *CI* (and ultimately on *DO*_{2}*i*), along with a Van’t Hoff dependence of *tVO*_{2}*i* on temperature (*T*):
Comparing the above to the equation of Figure 1, we see that, according to Equation 6, we must have *tVO*_{2} *i* = (1/0. 168) × 10^{0.0387T}, implying a Q10 of 10^{10×0.0387} = 2. 44.

### The GARIX+ Time-Series Model

The main assumption of the model is that changes (from one minute to the next) in the ‘logit’ of OER follow a homoscedastic (fixed-variance) normal distribution, and the distribution mean is a linear function of three ‘term groups’ (TGs): the disequilibrium TG (DTG), the autoregressive TG (ATG), and the exogenous TG (ETG). As a reminder, the logit transformation is given by:
It is useful for mapping a variable with a range between 0 and 1 (such as OER) to an unbounded range of − ∞ to + ∞, making it a more natural fit for linear models. The model, therefore, reads:
where the subscripts indicate measurement times (in minutes), *N* is the univariate normal (or Gaussian) distribution, μ_{t} is the time-dependent mean of the normal distribution, and σ^{2} is the time-independent variance of the distribution. The distribution mean at time *t*, μ_{t}, is given by:
The convention GARIX+(N) is used to refer to a GARIX+ model with a history length of *N* according to Equation 9. Also note that we have used the log transformation for state variables, *CI, Hb* and *SaO*2 in order to reduce their skewness.

Below we discuss each of the three TGs in the GARIX+ model. A conceptual summary of the model is provided in Figure 2.

### Disequilibrium Term Group (DTG)

The DTG represents a force to move OER in such a direction as to make VO_{2}i equal to tVO_{2}i. For this interpretation to be sensible, the parameter ‘k’ - which we refer to as the ‘disequilibrium coefficient’ - is expected to be positive: For instance, when VO_{2}i is larger than tVO_{2} i (over-oxygenation), then log(*VO*_{2} *i*_{t} ) − log(*tVO*_{2} *i*_{t} ) would be positive, and hence − *k*(log(*VO*_{2} *i*_{t} ) − log(*tVO*_{2} *i*_{t} )) would be negative (assuming *k >* 0). This induces a negative change in logit(OER), i.e., a decrease in OER, which would result in reduced VO_{2}i, thus bringing it closer to tVO_{2}i. The reverse would happen if VO_{2}i is smaller than tVO_{2}i (under-oxygenation).

Assuming that oxygen demand (tVO_{2}i) follows the Van’t Hoff parametric form, we have:
From the above, Q10 can be readily calculated:
Combining Equations 9 and 10, the following expression for DTG is obtained:
Since the above specification isn’t suitable for linear regression, we switch to the conventional specification:
which can be converted back to Equation 12 using:
Note that the above transformations are nonlinear, which means for estimating confidence intervals for *k, α*_{0}, *α*_{1}, we use Monte Carlo simulations, where samples are drawn from the joint normal distribution on (*η*_{0}, *η*_{1}, *η*_{2} ), returned by linear regression software in R, and *k, α*_{0}, *α*_{1} within each sample are computed according to Equations 14.

### Exogenous Term Group (ETG)

The ETG represents the effect of intended changes (the *n* = 0 terms in Equation 9) and past changes (the *n >* 0 terms in Equation 9) in the exogenous state variables (log(*CI*), log(*Hb*), log(*SaO*_{2} ), *T*) on change in *logit*(*OER*). Larger *n*’s correspond to more distant pasts. See Limitations for a discussion of endogeneity vs. exogeneity of state variables.

### Autoregressive Term Group (ATG)

The ATG represents the force that previous changes in OER exert on its future trajectory. Each term in this TG, β (*logit*(*OER*)_{n t−(n−1)} − *logit*(*OER*)_{t−n} ), corresponds to the change in logit(OER) from time *t* − *n* to time *t* − (*n* − 1), i.e., *n* minutes ago (since *t* represents ‘present’ time). If the coefficient, β_{n}, is negative, this term would produce a ‘self-correcting’ force, i.e., reversing stochastic changes in OER. For instance, if β_{1} < 0, then an increase in OER during *t* −*> t* + 1 would produce a force that seeks to reduce OER from *t* + 1 to *t* + 2.

When equilibrium is changed, for example due to cardiac index changing without any concomitant change in temperature, the force produced by the DTG is partially negated by the ATG contribution (Figure 3E). Note that ATG does not include an ‘intended’ term, i.e., *n* = 0, since this change corresponds to the outcome (endogenous) variable, which is not known at the time of making the prediction; rather, the change in logit(OER) from *t* to *t* + 1 is exactly what is being predicted.

### Static Model

As a reference point for the dynamic model, a naive, static model for learning the dependence of tVO_{2}i on temperature is considered. In this static model, an assumption is made that VO_{2}i is simply equal to tVO_{2}i at all times; in other words, both supply constraints and system dynamics are ignored. The relationship between tVO_{2}i and temperature is again assumed to follow the Van’t Hoff specification, as in the dynamic model. This leads to the following regression model:
which can be readily estimated from the perfusion data.

### Variable Importance Analysis

For Variable Importance (VI) analysis, a permutation-based method is used (Altmann et al, 2010), embedded in a repeated cross-validation (CV) scheme. In our context, the goal in conducting VI analysis is to quantify the relative importance of different variables, rather than variable selection. Permutation-based VI (PVI) is a model-agnostic VI method that measures degradation in predictive accuracy caused by random shuffling of the column (or columns) of the design matrix corresponding to the variables(s) in consideration. PVI is embedded in repeated cross-validation, in order to base the performance change on out-of-sample - rather than in-sample - predictive accuracy. Note that random shuffling applies not to the original time-series data but to the transformed version of the data that is used in estimating the linear regression coefficients. Relatedly, when applying PVI to TGs, the same random reordering of rows is applied to all individual terms in a TG (with each term corresponding to a column in the design matrix).

### System Simulations

Simulations are conducted to illustrate system response to exogenous changes. They can be done in two modes: deterministic and stochastic. In deterministic mode, the change in *logit*(*OER*) from time *t* to time *t* + 1 is taken to be the predicted mean, μ_{t}. In stochastic mode, the predicted mean is added to a random deviate drawn from a normal distribution with variance equal to σ^{2}, which is estimated during model training. Unless mentioned explicitly otherwise, initial conditions for simulations were based on median values of exogenous variables (*ci, Hb, SaO*2 and *T*) in the training data, and *OER* was calculated, using the GARIX+ model, such that the system would be in equilibrium (DTG = 0), where DTG is given by Equation 12.

## Results

### Patients and Outcomes

Table 1 shows a summary of patients used in the study and their outcomes. Note that, aside from gender and ethnicity, all other variables are per-operation. For instance, if a patient has had multiple operations, their age at the beginning of each operation is entered into the calculations. With regards to mortality, we see the observed 30-day mortality rate in our dataset is slightly below the expected mortality rate - for the same operations - calculated according to the PRAiS2 protocol.

### GARIX+ Model Fit

Figure 3 shows the results of fitting the above-described GARIX+ model to the granular perfusion data. The magnitude of coefficients for autoregressive and exogenous terms declines with lag increase and reaches statistically insignificant levels after 10 minutes or so, consistent with a finite memory in the system (Panel A). (See Supp Mat E for the complete list of coefficients for the GARIX+(20) model.) Coefficient signs are consistent with an adaptive OER change to: 1) maintain oxygen consumption at a fixed temperature (negative coefficients for change terms involving log(*CI*) and log(*Hb*)); 2) increase/decrease oxygen consumption when temperature is increased/decreased, respectively (positive coefficients for change terms involving *T*); and 3) stabilise random fluctuations in oxygen consumption (negative coefficients for *logit*(*OER*) change terms). *SaO*_{2} is a special case since it is also part of the OER definition: *OER* = 1 − *SvO*_{2}/*SaO*_{2}, and it therefore influences VO_{2}i both directly, and also via OER:
As shown in Methods, the DTG coefficients can be used to calculate Q10 and construct the Van’t Hoff curve (Panel B). The estimated Q10 is 2.2, which is close to the 2.4 figure given in Kirklin & Blackstone (2012), and significantly below the 3.7 value estimated from the data using a naive, i.e. static, model (see Static Model). Also, ‘k’, the disequilibrium coefficient, is estimated to be 0.024, with a 95% confidence interval of 0.021-0.026 that is significantly positive. Positivity of the disequilibrium coefficient is consistent with an adaptive force to bring VO_{2}i in line with target VO_{2}i.

By construction, the static fit assumes that the data is balanced between episodes of under- and over-oxygenation. (Note, however, that the static regression is performed on the logarithm of VO_{2}i, which means residuals may not be balanced on the original scale, especially if the model is misspecified.) Yet, within each subject, and especially at the granular, minute-by-minute level, there are significant deviations from the regression line, as can be seen in Figure 3B. On the other hand, the GARIX+(7) fit implies more cases of over-oxygenation, especially at higher temperatures. While the ‘Kirklin - asymptote’ line has a Q10 close to GARIX+, the absolute levels are significantly higher. The discrepancy becomes smaller when we use the ‘clinical’ values given in Figure 2-11 of Kirklin and Blackstone (2012). Either way, the Kirklin lines imply under-oxygenation to be a more prominent risk, compared to the fit obtained from our GARIX+ model. Panel C confirms that only 5-10 minutes of history are useful for improving the out-of-sample predictive accuracy of the model. This observation is the basis for using GARIX+(7) in most of the remaining analyses in this paper. It is possible, however, that larger training data would justify using a longer history.

Panel D shows PVI by TGs (DTG, ATG, and ETG). Note that DTG - which allows the estimation of tVO_{2}i and its dependence on temperature - has the smallest, but statistically significant, contribution of the three TGs. This limits the ability to conduct subgroup analyses, and to include many patient attributes and interaction terms in the DTG of the model.

Panel E shows the timeline of contributions from different TGs to change in *logit*(*OER*), driven by a hypothetical step-function change in CI at t = 0. The initial, fast adaptation is dominated by the ETG, while the secondary, slow adaptation is dominated by the counteracting effects of the DTG and the ATG. The ATG contribution (red) is partially offsetting the DTG (green), resulting in a slower overall adaptation rate (purple). Panel F illustrates the resulting OER trajectory consisting of the initial fast stage, followed by the secondary, slow stage.

### Subgroup Analysis by KDIGO Score

Next, we trained the GARIX+(7) model on subsets of patients that all share a single KDIGO score. Comparing these models allows us to identify interesting correlates of AKI severity.

In Panel A of Figure 4, the model fit, measured by the out-of-sample (OOS) R-squared, degrades with increased AKI severity (i.e., increasing KDIGO score; red line labelled ‘Genuine’), and this decline cannot be explained by changes in the sample size (blue line, labelled ‘Random’). This is directionally consistent with the model reflecting, at least to some degree, normal body mechanisms (which would presumably be stronger in AKI-free patients), rather than a purely pathological set of interactions triggered by the extreme conditions that patients undergo during CPB.

Panel B shows that the immediate adaptive reaction to changes in haemoglobin is an important factor in explaining the change in model performance with AKI: Those with severe AKI show a much less pronounced immediate adaptation to a change in blood haemoglobin levels. Also, as panel C shows, the disequilibrium force becomes steadily weaker (smaller disequilibrium coefficient or ‘k’) in patients with more severe AKI, which is, in turn, reflected in a slower adaptation time (Panel D). Finally, smaller Q10 values in higher-KDIGO patients (corresponding to the *α*_{1} parameter in Panel C) suggests a weaker dependence of metabolic needs - genuine or pathological - on body temperature in severe-AKI cases.

### Comparison to the Kirklin Model

The elbow, or piecewise linear, model - where VO_{2}i becomes completely supply-independent after DO_{2}i exceeds a threshold - is appealing due to its simplicity, and also because it provides exact estimates not only for tVO_{2}i, but also for tDO_{2}i.

However, empirical evidence (Kirklin & Blackstone, 2012) suggests that a hyperbolic model of oxygen consumption is more accurate in explaining the - albeit limited and dated - human and animal data.

The GARIX+ dynamic model offers an alternative explanation of empirical data. Figure 5 shows the results of simulating the human experiments reported by Fox et al (1982).

With a wait time of 10 minutes between state transitions, GARIX+(7) model produces a DO_{2}i-VO_{2}i curve that is similar to the hyperbolic fit of Kirklin. As wait time is increased, the GARIX+(7) model produces sharper curves, and at infinity it converges to the piecewise linear shape. As shown earlier (Figure 3B), the predicted tVO_{2}i is higher in Kirklin than GARIX+(7) at all temperatures. As shown in Table 2, the distributions of R-squared, ‘a’ and ‘b’ parameters (using Fox et al convention, see Supp Mat B) according to GARIX+(7) - and using a wait time of 10 minutes between state transitions - all have 95% confidence intervals that cover the results reported in Fox et al.

While a static view of the piecewise linear model would only allow for the possibility of under-oxygenation, the shift to a dynamic perspective opens up both possibilities. This is illustrated in Figure 6, which shows a 3D surface plot of oxygenation gap (OG) (y axis) vs. the relative step change in DO_{2}i at t = 0 (x axis), and at different times since the introduction of step change (z axis). OG is defined as the % difference between oxygen consumption at time *t* (*VO*_{2} *i*_{2} ), and target oxygen consumption (*tVO*_{2} *i*):

*OG*_{t} = 100 × (*VO*_{2} *i*_{t} − *tVO*_{2} *i*) / *tVO*_{2} *i*. (Note that tVO i does not change with time since temperature is fixed.) It can be demonstrated that, for small *t*’s, the response is nearly linear, corresponding to relatively little adaptation in OER. This means that, e.g., when cardiac index is increased and before OER adaptively decreases to bring VO_{2}i down to its pre-perturbation level, there will be a nearly proportional oxygen consumption surplus, or over-oxygenation. However, as time progresses, such oxygen surplus dissipates as OER adaptively decreases. While oxygen surplus eventually vanishes, oxygen deficit cannot completely disappear if the ratio *tVO*2*i*/*DO*2*i* implies an OER greater than 1. As a result, the piecewise linear shape emerges for very large *t*’s.

### Q10 by Temperature Band

One of the limitations of the Kirklin model (Equation 6) is that the dependence of tVO_{2}i on temperature was assumed to follow the strict parametric form of Van’t Hoff (where the logarithm of oxygen demand is a linear function of temperature), and yet the evidence for this strong parametric assumption is scant. In fact, Equation 6 is based only on two data points: a set of human experiments conducted in 1982 at 20°C (Fox et al, 1982), and three sets of animal experiments conducted in the 1950s, at 37°C (Cheng et al, 1959; Paneth et al, 1957; Starr, 1959). While the Van’t Hoff specification has been followed in Equation 9 due to its widespread adoption, it may be informative to perform a subgroup analysis by temperature in order to see if the GARIX+(7) model arrives at the same, or close values of Q10 within different temperature bands.

As Figure 7A shows, most of the temperature recordings are clustered in the 30-37°C range. This, combined with a low relative VI for the DTG means that models trained within fine temperature subgroups may not yield statistically significant results for the parameters of DTG. Therefore, we opted to create two temperature bands, separated at 28°C (Panel B). We observe that, GARIX+(7) models that are fit within each subgroup yield significantly different oxygen demand curves. While the low-temp regime has a Q10 of 2.6, the high-temp regime produces a Q10 that is barely above 1, suggesting a much less pronounced dependence of oxygen demand on temperature in this range.

## Discussion

This study is the first to present a data-driven, analytical framework for predicting OER changes in response to clinical interventions during CPB. The analysis was made possible through a combination of granular perfusion data, combined with two important methodological enhancements to an otherwise standard framework for time-series analysis.

Our model suggests that changes in the components of oxygen delivery cause an adaptive OER response to keep oxygen consumption in line with oxygen demand, according to body temperature. The dependence of oxygen demand on temperature decreases as temperature increases towards normothermia, calling into question the accepted Van’t Hoff equation. Children developing AKI exhibit a dampened response to changes in haemoglobin during CPB, with this dampening of response intensifying with AKI severity. This association can be the foundation of an intraoperative tool for early diagnosis of at-risk patients. Furthermore, our model provides a novel framework for GDP by identifying periods of over- and under-oxygenation through comparing the VO_{2}i against tVO_{2}i, and by predicting the impact of perfusionist actions on OER and oxygen consumption.

### GARIX+ as a Grey-Box Model

In creating GARIX+, a standard autoregressive, integrated time-series model with exogenous variables (ARIX) - a black-box model - has been enhanced in two important ways. First, a DTG has been added (hence the ‘+’ sign) that reflects our expectation that the body seeks to keep oxygen consumption in line with its metabolic needs, as dictated by body temperature. Secondly, the data points are combined across all patients, resulting in a ‘global’ time-series model. This gave the model sufficient statistical power to estimate the DTG parameters, despite their relatively-small VI compared to the ATG and ETG. Within the DTG, the dependence of oxygen demand (tVO_{2}i) on temperature is assumed to follow the Van’t Hoff model, which is the standard specification in the literature. As such, GARIX+ can be characterised as a ‘grey-box’ model; that is, a black-box model enhanced with domain knowledge.

In a GARIX model (no ‘+’), any combination of state variables (CI, Hb, SaO_{2}, OER and temp) would be acceptable as a solution, as long as they all remain constant with time. While a GARIX model could explain much of the OER variations in the data (as evidenced in low relative VI of the disequilibrium TG in Figure 3D), from a physiological perspective it would be an incomplete model since it does not reflect the expectation that the system would seek to adapt oxygen consumption to match the metabolic needs of the body. At the same time, the linear specification of GARIX+, which was our primary reason for not using machine learning models, allowed us to examine and interpret the model coefficients and gain confidence in the physiological plausibility of the results, a topic that is further discussed next.

## Physiological Interpretation of Findings

### The GARIX+ Model

Oxygen, whilst vital, can also be toxic, with both oxidative stress and cellular hyperoxia identified as major contributors to multiple pathological states (Evans, 2016). Historically, under-oxygenation was (and still remains) the key concern in patients undergoing cardiac surgery, although recent studies are starting to highlight the pathological dangers of over-oxygenation during CPB, with an ongoing debate around the level of damage (Jakutis et al, 2017; Young 2012; Abou-Arab et al, 2020). Equally, being able to define what these terms (over- and under-oxygenation) mean is complex. For example, within the kidney’s microcirculation, the renal medulla is normally exposed to lower partial pressures of oxygen compared to the renal cortex, therefore under- and over-oxygenation will be very different conditions, even in the same organ (Evans et al, 2008). Hyperoxia, however, is a phenomenon not encountered in the natural world, and therefore one could speculate that there has been little evolutionary need to develop countermeasures to it (although atmospheric oxygen levels have historically been much higher than they are today; Lane 2003). At the cellular level, oxygen is involved in the final step (Electron Transport Chain, ETC) of aerobic metabolism, in the mitochondria. Whilst the citric acid (Kreb’s) cycle will only occur in the presence of oxygen, it will continue even to a mitochondrial PO_{2} level as low as 3mmHg (0.4 kPa), with a *K*_{M} of 22 mmHg (1.6 kPa; Wilson et al, 2012).

Therefore, an explanation for the piecewise linear model is that once sufficient oxygen molecules are available (oxygen consumption) to support the ETC, further increases in partial pressure will not increase the rate of metabolism. Moreover, as oxygen concentrations reach excessive levels, animal and cell culture models have demonstrated diminished ETC function and altered mitochondrial morphology (Resseguie et al, 2015). Paradoxically, hypoxia is also associated with an increased release of Reactive Oxygen Species (ROS) and mitochondria themselves appear to be O_{2} sensing although our understanding of the mechanisms behind this is still incomplete (Guzy and Schumacker, 2006). Hence, there are clear indications that a balance in oxygen consumption and extraction is playing out at the cellular level, and supports our notion for a Disequilibrium Term Group at the systemic level; i.e., the body is attempting to maintain VO_{2} in line with tVO_{2}. Given that individual tissues and organs are known to operate at varying PaO_{2} levels or with different metabolic requirements, and even require particular OER levels to maintain proper concentration gradient of oxygen, it follows that, at a systemic level, alterations to exogenous variables (represented in the Exogenous Term Group in the model) will have shorter or longer effects in those tissues (Wolff, 2013). This lag or memory of previous events may be also partially a result of perfusionist interventions coming in relatively rapid succession and independently of normal physiological mechanisms. Such deviations from the normal autoregulatory patterns are frequently seen during CPB, and indeed, a recent review indicated that cardiac surgery is associated with changes in cerebral autoregulation which significantly relate to clinical outcomes (Caldas et al, 2018).

### Dependence of Q10 on Temperature Regime

Oxygen can regulate glycolysis via the Pasteur effect. The availability of oxygen will suppress glycolysis, and decreased availability leads to an acceleration of glycolysis, at least initially. This appears to be regulated by several factors including the allosteric regulators of glycolysis (enzymes such as hexokinase; Lenzen 2014).

Allosteric enzymes do not follow the Van’t Hoff dependence on temperature, instead often displaying a sigmoid shape for dependence of reaction velocity on temperature, at a fixed substrate concentration.

Whilst Q10 has been traditionally used to describe the dependence of rate of chemical reactions on temperature, enzymatic reactions have ranges over which the proteins can operate before conformational changes or denaturing occurs. Furthermore, bypass physiology is very different from normal human physiology, in that exogenous variables (such as haemoglobin concentration, cardiac index etc.) are manipulated separately by the perfusionist, rather than in coordination, as would be under normal biology. For instance, normal physiological mechanisms to combat hypothermia (such as shivering) are prevented through general anaesthesia, and PaO_{2} is managed at levels higher than observed naturally, in order to allow a margin of safety. In the analysis described in this study, a lower Q10 value was measured in patients that were treated with mild hypothermia (>28°C), compared to those with moderate and deep hypothermia (<28°C). One might postulate that under mild hypothermic conditions sufficient substrate and oxygen is being provided, and enzymatic functions are maintained such that they lead to a degree of decoupling between metabolic rate and temperature. On cooling further into moderate and deep hypothermia, this independence is lost and we thus return to the more traditional Q10 of 2-3. However, further, independent data are required to validate these observations and experimentation is needed to determine the cellular mechanisms operating in different temperature regimes.

### Altered Model Performance in Patients Developing AKI

The aetiology of AKI is complex and incompletely understood. Some studies have pointed to renal haemodynamics and oxygenation as AKI triggers (Evans et al., 2018), whilst a meta-analysis suggested hypertension, alterations of serum creatinine, CPB and aortic clamping times (among others), as strong predictors for the development of postoperative AKI (Yi et al 2016). Renal oxygenation is directly dependent on renal blood flow and its regulation, and its oxygen consumption levels are primarily driven by tubular demand, which itself is dependent on the glomerular filtration rate (Smith et al 1940). There exist, within the kidney, regional variations in oxygen delivery and consumption; areas with normally lower physiological pO_{2} are at greater risk of hypoxic damage through alterations in the delivery - consumption balance (Edwards and Kurtcuoglu, 2022).

The reduced performance of the GARIX+ model in severe AKI patients, as well as their slower adaptation seems physiologically plausible, as it suggests a disruption of feedback pathways in severe-AKI patients. The causality mechanism is not clear, however. For instance, slower adaptation may be causing AKI by creating prolonged periods of imbalance between oxygen consumption and demand, or there may be a common cause for slower OER adaptation during CPB, and postoperative AKI. One of the limitations of this (and indeed nearly every other) study examining AKI after CPB, is a lack of direct measurement of blood flow and oxygenation occurring in the kidney. Therefore it is impossible to provide a definitive explanation for the slowing of the adaptive response in high-severity AKI patients (Figure 4D). However, this is consistent with an exaggerated supply-dependence of consumption in critically-ill patients, including those with sepsis (Friedman et al, 1998).

An intriguing observation from the GARIX+ model, is the decline in adaptive response to haemoglobin changes in severe-AKI patients. Furthermore, although not detailed in this manuscript, a progressively weaker haemoglobin response was observed in younger patients. During CPB, other than through haemofiltration, the only way to increase Hb concentration is through the addition of stored RBCs. Due to the large volume of the CPB circuit relative to the child’s circulating volume, RBC transfusions are more common in paediatric than adult surgery. A recent study examined the effect of RBC transfusion on renal oxygenation in kidneys supported via an *ex vivo* normothermic kidney perfusion rig (Bumbill et al 2023). It demonstrated that there was slower O_{2} unloading from stored RBC with concomitant reduction in O_{2} extraction in perfused kidneys, indicating diffusion-limited O_{2} release at capillaries. It is also possible that SaO_{2} dependence on PaO_{2} - i.e. the oxygen dissociation curve - shifts differently with transfused RBC’s due to the metabolic and morphological changes that can occur due to storage (Bennett-Guerrero, 2007 & Yoshida, 2019). Further examination is required to ascertain the underlying mechanism for this response, but it presents a potential mechanism for early - i.e., intraoperative - diagnosis of patients at risk of post-surgical AKI.

### Implications for Goal-Directed Perfusion

Whilst the majority of research has focused on preventing under-oxygenation of the tissues during CPB, little has been done to date to assess over-oxygenation. Furthermore, DO_{2}i provides limited information and assumes an optimised dissociation curve occurring across all tissues. Using predicted tVO_{2}i as a reference (see Table 3) will allow the more specific definition of over-and under-oxygenation, increasing the predictive power towards AKI and other adverse events.

## Limitations

### Contribution from dissolved oxygen

The contribution from dissolved oxygen was not considered due to a lack of access to PvO_{2} (partial pressure of mixed venous oxygen). However, as stated before, this contribution is on average only around 1.5% of the total oxygen content. Therefore, the analysis and results can be considered directionally correct. Future work will incorporate the dissolved-oxygen component and improve the accuracy of the models.

### Endogeneity vs. Exogeneity of State Variables

Exogeneity of CI, Hb, SaO_{2} is an idealised assumption. In reality, values of these state variables at time “t+1” cannot be set with full certainty at time “t”; rather, these variables are controlled indirectly. In other words, the boundary between exogeneity and endogeneity is blurred. In particular, SaO_{2} is controlled via PaO_{2} (changing the arterial oxygen pressure). A basic dynamic model confirmed that PaO_{2} (plus temperature) can predict changes in SaO_{2} with good accuracy (results not shown here). In future work, more realistic models will be needed to better reflect the true controllability of the ‘exogenous’ variables.

### Constraints on OER

Other than using a logit transformation to implicitly constrain OER to fall between 0 and 1, there were no other constraints on OER. The nonlinearity of the logit transformation also provides a natural slowing down of OER change as it approaches extreme values, which lends further physiological plausibility to the system simulations. It is possible that, during normal conditions, the body deploys other adaptive mechanisms, which may or may not be disrupted during CPB, to ensure that OER remains within a tighter range. During CPB, on the other hand, perfusionists generally go to great lengths to ensure that OER does not stray too far from an acceptable range. For instance, in our data, 90% of time OER remains below 31%. As such, during system simulations, care must be taken while interpreting results where simulated OER approaches or exceeds the observed range in our data.

### Paediatric patients

This dataset contains a broad paediatric age range and yet no adults. Mixing neonates and teenagers, who may have very different physiology and different responses in terms of their adaptation to CPB stresses, is not ideal. However, due to the relatively low VI for the DTG (OER variations are dominated by exogenous and autoregressive forces during CPB), we decided to combine all patients to build a single, global model in order to maximise the study power with regards to the steady-state analysis (dependence of tDO_{2}i on temperature). More data and analyses are needed to differentiate between age groups.

## Conclusion

This study proposes a physiologically plausible model of OER changes during CPB that is consistent with past experimental data. While in the long-term under-oxygenation is the dominant risk, slow adaptation of OER during CPB creates short-term opportunities for over-oxygenation after significant changes in the exogenous variables, such as a sudden increase in cardiac index. The model provides well-defined values for tVO_{2}i at a given temperature (Table 3), paving the way for further research into the effects of over- and under-oxygenation on postoperative outcomes such as AKI, and hence improvements in GDP protocols. Furthermore, the identification of an altered response to haemoglobin changes in patients that proceed to severe post-surgical AKI suggests a potential intraoperative diagnostic tool for early identification of at-risk patients.

## Data Availability

All data produced in the present study are available upon reasonable request to the authors

## Supplementary Material

## A: Approximate Equivalence of Van’t Hoff and Arrhenius Specifications

The Arrhenius and Van’t Hoff formulations are the best known parametric forms used to describe the relationship between temperature and oxygen consumption. It can be demonstrated that, as long as temperature changes are small compared to their absolute values, *when expressed in Kelvins*, the two specifications are approximately equivalent. In this manuscript, the Van’t Hoff specification has been adopted, which is consistent with the widely recognised concept of Q10, the multiplicative increase in RMR for every 10 degrees C increase in body temperature.

The Arrhenius equation, after mapping the concept of chemical reaction rate to metabolism and hence oxygen consumption, would require the logarithm of oxygen consumption to be a linear function of the inverse of temperature, expressed in Kelvins ( *T*_{k} ):
Let’s express *T* in terms of its deviations (*ΔT*) from a baseline value, *T*_{0} :
Equation (A.1) can be rewritten as:
If *ΔT* << *T*_{0}, i.e., if changes in temperature are small compared to its baseline value, in

Kelvins, then a first-order, Taylor-series expansion is used to rewrite Equation (A.3):
Note that the last expression on the right is a linear function of temperature change, *ΔT*. This is simply the Van’t Hoff specification, i.e., logarithm of VO_{2}i being a linear function of temperature. (Keep in mind that one Kelvin change equals one Celsius change in temperature, i.e., *ΔT* would be identical whether expressed in Kelvins or Celsius.)

## B: Experimental Setup and Model Specification of Fox et al (1982)

Table 1 in Fox et al (1982) was consulted for constructing the random sequence of CI values for each simulated subject who received 4 step-function changes to their CI during the simulations. (The initial value of CI was 2.2 L/min/m^{2} for all subjects.) This means that, for each step-function change in CI, one of the four possible CI bands listed in Table 1 is first randomly selected and then a random number is drawn from the normal distribution associated with that CI band, according to the mean and standard deviation values listed therein. Fox *et al*., fitted a hyperbolic model to their experimental data, which assumes 1/VO i to be a linear function of 1/CI: *VO*_{2} *i* = *a* × *CI* / (*b* + *CI*), with *a* = 35, *b* = 0. 42. Note that this equation is equivalent to Equation 6 for *T* = 20*C*.

## C: GARIX+ Model Diagnostics

Panel A of Figure C-1 shows the combined autocorrelation function (ACF) plot of GARIX+(7) model residuals across all segments. For comparison, Panel B shows the same plot, where model residuals are replaced with white noise. The overall similarity of plots - with the average (red) line being close to the y = 0 line - indicates that model residuals are well-behaved. Panel C shows the ACF plot for the absolute of model residuals, and Panel D shows the ACF of the absolute of white noise. The discrepancy between these two panels is a sign that there is information in the magnitude of residual noise, i.e., heteroscedasticity. Incorporating heteroscedasticity into the GARIX+ model is a topic of future research.

## D: Statistical Significance of GARIX+ Coefficients

The supplementary file, GARIX_20_model_coefficients_supp_mat.csv, shows the list of coefficients and their statistical significance for the GARIX+(20) model. The majority of terms have statistical significance for up to 5-10 minutes of lag, and the terms corresponding to recent history have particularly high significance. An exception is SaO_{2}, which is only highly significant at lag = 0 (corresponding to change from t to t+1).

## Footnotes

↵

^{1}Note that the hyperbolic function is mathematically equivalent to the Michaelis-Menten formulation, a terminology used for describing enzymatic reactions.

## Nonstandard Abbreviations and Acronyms

- AKI
- acute kidney injury
- GARIX
- global autoregressive Integrated (time-series model) with exogenous variables
- GARIX+
- GARIX model, plus a disequilibrium term group (in addition to autoregressive and exogenous term groups)
- ATG
- autoregressive term group
- BSA
- body surface area (m
^{2}) - CI
- cardiac index (blood flow indexed BSA) (L/min/m
^{2}) - CPB
- cardiopulmonary bypass
- CV
- cross-validation
- DHCA
- deep hypothermic circulatory arrest
- DO
_{2} - oxygen delivery (mL/min)
- DO
_{2}i - oxygen delivery, indexed to BSA (mL/min/m
^{2}) - DTG
- disequilibrium term group
- EHR
- electronic health record
- ETG
- exogenous term group
- GDP
- goal-directed perfusion
- Hb
- haemoglobin content of blood (g/dL)
- HLM
- heart-lung machine
- i.i.d.
- Independent and identically distributed
- ML
- machine learning
- OER
- oxygen extraction ratio
- OG
- oxygenation gGap
- OOS
- out-of-sample
- PaO
_{2} - partial pressure of oxygen (kPa)
- PVI
- permutation-based variable importance
- Q10
- multiplicative increase in the resting metabolic rate for every 10°C increase in body temperature
- RBCs
- red blood cells
- RMR
- resting metabolic rate
- RTO
- return-to-equilibrium (an estimated parameter in the GARIX+ model)
- SaO
_{2} - oxygen saturation of arterial blood (scale: 0-100)
- SvO
_{2} - Oxygen saturation of mixed venous blood (scale: 0-100)
- TG
- term group
- tVO
_{2}i - target (indexed) oxygen consumption, aka (indexed) oxygen demand (mL/min/m
^{2}) - tDO
_{2}i - threshold (indexed) oxygen delivery (mL/min/m
^{2}) - VI
- variable importance
- VO
_{2} - oxygen consumption (mL/min)
- VO
_{2}i - oxygen consumption, indexed to BSA (mL/min/m
^{2})