Temporal Changes in the Oxyhemoglobin Dissociation Curve of Critically Ill COVID-19 Patients

Critical COVID-19 is a life-threatening disease characterized by severe hypoxemia with complex pathophysiological mechanisms that are not yet completely understood. A pathological shift in the oxyhemoglobin curve (ODC) was previously described through the analysis of p50, intended as the oxygen tension at which hemoglobin is saturated by oxygen at 50%. The aim of this study was to analyze Hb-O2 affinity features over time in a cohort of critically ill COVID-19 patients, through the analysis of ODC p50 behavior. A retrospective analysis was performed; through multiple arterial blood gas (ABG) analyses, each p50 was calculated and normalized according to PaCO2, pH and temperature; patients’ p50 evolution over time was reported, comparing the first 3 days (early p50s) with the last 3 days (late p50s) of ICU stay. A total of 3514 ABG analyses of 32 consecutive patients were analyzed. The majority of patients presented a left shift over time (p = 0.03). A difference between early p50s and late p50s was found (20.63 ± 2.1 vs. 18.68 ± 3.3 mmHg, p = 0.03); median p50 of deceased patients showed more right shifts than those of alive patients (24.1 vs. 18.45 mmHg, p = 0.01). One-way ANOVA revealed a p50 variance greater in the early p50s (σ2 = 8.6) than in the late p50s (σ2 = 3.84), associated with a reduction over time (p < 0.001). Comparing the Hb-O2 affinity in critically ill COVID-19 patients between ICU admission and ICU discharge, a temporal shift in the ODC was observed.


Introduction
The oxyhemoglobin dissociation curve (ODC) is a sigmoid curve representing the relationship between arterial hemoglobin saturations (SaO 2 ) at different arterial blood oxygen tensions (PaO 2 ) [1,2]; the continuous correlation between SaO 2 and PaO 2 is reported as hemoglobin-oxygen (Hb-O 2 ) affinity. ODC was first described by Christian Bohr in 1904 [3], who demonstrated how this affinity was mathematically represented by a sigmoid function.
p50 is the oxygen tension at which hemoglobin is saturated at 50%; in normal conditions, it ranges from 24 to 28 mmHg [4]. It follows the fact that a specific p50 value corresponds to a specific ODC curve and, therefore, to a specific Hb-O 2 affinity. Various biochemical variables can physiologically influence the Hb-O 2 affinity, leading to an ODC shift and to a consequent p50 change [5], with right and left shifts occurring in relation to changes in blood pH, PaCO 2 , temperature and 2,3-biphosphoglycerate (2,3-BPG) plasma concentration (Bohr/Haldane effect). The ODC shift can also be induced by other physiological or pathological conditions, such as aging [6], diabetes mellitus [7], anesthesia [8] and respiratory failure [9]. Considering these variables, an equation to represent the corrected ODC and its p50 was developed. In 1910, Archibald Vivian Hill, a British physiologist, was the first to elaborate on this equation [10], which was then revised by John Severinghaus in 1979 [11] and subsequently by many authors. Dash et al. [12] proposed a simple accurate mathematical model based on the original Hill formula with variation for O 2 , CO 2 , pH, 2,3-BPG and temperature ( Figure S1).
COVID-19 is a disease characterized by severe hypoxemia [13] with complex pathophysiological mechanisms that are not yet completely understood [14]. In the acute phase, tachypnea and hypoxemia usually develop; quite surprisingly, and in contrast with similar conditions affecting the respiratory physiology, severe hypoxemia in COVID-19 is not always associated with dyspnea [14][15][16][17]. Some authors hypothesized that this discrepancy may be due to ODC shifting in the context of SARS-CoV-2 infection, possibly due to direct viral modulation on Hb structure, as evidenced by Wenzhong et al. [18] and Lanxioux et al. [19]. Vogel et al. [20] performed a retrospective analysis on arterial blood gas (ABG) analyses of more than 3500 blood gas analyses in 43 critically ill COVID-19 patients, comparing them with historical controls of patients with acute hypoxemia but without SARS-CoV-2 infection, observing a left shift of the ODC, whereas other authors did not find any changes [21,22]. Their results therefore suggest an in vivo increase in Hb-O 2 affinity in critically ill COVID-19 patients. In contrast, Daniel et al. [21] observed that the in vitro Hb-O 2 affinity in 14 patients affected by COVID-19 did not differ from that of 11 control patients in standard conditions; similarly, Pascual-Guardia et al. [22] did not find any change in the Hb-O 2 affinity in COVID-19 patients.
Taking into consideration the temporal evolution of COVID-19 [23,24], and the feasible modification of the Hb structure induced by SARS-CoV-2 itself [18,19], we hypothesized temporal changes in the ODC shift in critically ill COVID-19 patients, which possibly correlates with pulmonary gas exchange, disease progression and survival. The aim of the study was to analyze Hb-O 2 affinity features and their modifications over time through the analysis of normalized p50 behavior in a cohort of critically ill COVID-19 patients, hoping to aid in unveiling the underlying pathophysiological mechanism of this emerging disease in a pandemic scenario.

Materials and Methods
This was a retrospective observational study in a cohort of consecutive COVID-19 patients admitted from March to May 2020 to our COVID-19 Center in Lugano, Switzerland. The inclusion criteria were critically ill COVID-19 patients older than 18 years of age of both sexes undergoing invasive mechanical ventilation; diagnosis of COVID-19 was made by a positive RT-PCR molecular test. Patients excluded from the study were those with less than seven ABG analyses in the health record or with less than seven days in ICU care. Demographics and clinical and laboratory data were obtained from the hospital's health medical record.
Each ABG analysis was carried out using one of the four ABL90 FLEX blood gas analyzers (Radiometer©) available in the intensive care unit (ICU) and emergency department (ED). Parameters such as the pH, PaO 2 , PaCO 2 , SaO 2 , sodium (Na + ), potassium (K + ), calcium (Ca 2+ ), chloride (Cl − ), glucose, L-lactate and hemoglobin were directly measured by the machine. Temperature and FiO 2 at the time of arterial blood gas analysis were instead manually reported from the blood gas analyzer by the nurse in charge. Through specific queries already available in the health medical records, all patients' ABG analyses, performed during the study period, were obtained, each of them reporting patients' identification number and date and time of the analysis, as well as the pH, PaCO 2 , PaO 2 , SaO 2 , Na + , K + , Ca 2+ , Cl − , glucose, L-lactate and hemoglobin values. All the above-mentioned data were collected in an Excel dataset.

p50 Calculation with Variable Coefficient
We used the Dash formula [12] to calculate the ODC p50, starting from the PaO 2 /SaO 2 relationship through the modified Hill coefficient according to Dash [10,12,25] and setting the parameters α, β and γ to 2.82, 1.2 and 29.25, respectively, in accordance with Severinghaus' experimentally available data [11]. A weighted correction, balanced for PaCO 2 , pH and temperature, was implemented using the same Dash formula and the same definition as above, considering 2,3-BPG at 4.65 mmol/L as a constant ( Figure S1) [1,4,20,[26][27][28]. It was then possible to obtain a weighted p50 for each corresponding arterial blood gas analysis' PaO 2 /SaO 2 pairing, allowing a comprehensive analysis and comparison of the derived oxygen dissociation curves. According to the mathematical limitation in the Hill function [12], consisting of an excessive p50 increase when PaO 2 reached high values (above 100 mmHg), we excluded all arterial blood gas analyses with PaO 2 greater than 100 mmHg ( Figure S2).
First, the standard ODC was calculated according to the original computations performed by Severinghaus et al. [11]. Based on the various modifications of patients' p50 values occurring during the ICU stay, which were weighted according to the Hill and Dash formula [12], each corresponding ODC was then graphically reported and compared to the original one.

The p50 Shift over Time
For each patient, we described the p50 shift over time, plotting p50 values for each arterial blood gas analysis and thus representing the ODC temporal change. Moreover, to better characterize how relevant differences correlated with the disease severity evolution during the ICU stay, we calculated each patient's p50s from the first three days in the ICU (early p50s) and compared them with those from the last three days in the ICU (late p50s). The time frame of three days was arbitrarily chosen as a balance between obtaining a high number of blood gas analyses and, at the same time, avoiding analyzing too many days of treatment altogether, especially given that numerous factors could rapidly alter the patients' clinical condition (bacterial superinfection, sepsis, etc.).
Patients were further stratified according to their own temporal shift in "right-shift p50" and "left-shift p50", and a further subgroup analysis was performed to identify any specific correlation, especially regarding p50 change over time. An additional comparison between all p50 in alive and deceased patients was performed. Finally, starting from the preliminary results, a large variability in the p50 of critically ill COVID-19 patients over time was identified. To better characterize this aspect and to define whether the p50 temporal evolution is constant and regular or fluctuates over time, a variance analysis for each single patient during all ICU hospitalizations was performed, which further analyzed any eventual variance in temporal progression.

Outcomes
The primary outcome was to assess the temporal p50 change in critically ill COVID-19 patients, comparing early p50s vs. late p50s, to describe the Hb-O 2 affinity change over time. Secondary outcomes were to correlate the ODC and p50 variabilities with ICU survival and to stratify patients according to the temporal patterns of different p50 shifts, to analyze the p50 variability over time.

Statistical Analysis
Descriptive statistics were calculated; categorical data were reported as numbers (percentages). Numerical data distribution was reported as the mean (SD) or as median (IQR) according to statistical distribution, verified by the Kolmogorov-Smirnov test. Differences between continuous data were studied with the t-test for independent groups or with the Wilcoxon test if nonparametric analysis was required. The study of differences between the groups for categorical data was carried out with the chi-square test. A variance analysis was carried out with a one-way ANOVA through the clinical variables-sex, age, BMI, ICU length of stay (LOS), arterial hypertension, diabetes, chronic obstructive pulmonary disease (COPD) and obstructive sleep apnea syndrome (OSAS)-the prognostic variables--Simplified Acute Physiology Score (SAPS II), Sequential Organ Failure Assessment (SOFA) and the nine equivalents of nursing manpower use score (NEMS scores)-and laboratory variables-lactates, leucocytes, lymphocytes, thrombocytes, aspartate aminotransferase (AST), alanine aminotransferase (ALT), total bilirubin, C-Reactive Protein (CRP), lactate dehydrogenase (LDH), creatine kinase (CK), ferritin and creatinine; one-way multivariate analysis of variance (MANOVA) was also performed. All differences were reported as 2-tailed; an equal variance was not assumed. The confidence interval (CI) was 95%; the level of significance was established to be <0.05. No imputation for missing data was carried out. Statistical data analysis was performed using the SPSS 26.0 package (SPSS Inc., Chicago, IL, USA).

Results
During the enrollment period, 40 consecutive patients were admitted to the ICU, and 3747 total blood gas analyses were performed; 3 (7.5%) patients were excluded due to an ICU LOS less than 7 days, 5 (12.5%) patients were excluded due to a number of ABG analyses less than 7, 210 ABG analyses were excluded for a PaO 2 greater than 100 mmHg. A total of 3514 ABG analyses from 32 patients were therefore analyzed ( Figure 1). A total of 29 (87.8%) patients were men, with a mean age of 62 years (SD 11.7); 12 (36.3%) patients were affected with diabetes, 14 (42.4%) with arterial hypertension, 3 (9%) with COPD, 5 (15.1%) with OSAS and 5 patients (15.1%) presented with a diagnosis of pulmonary embolism at ICU admission. The mean SAPS II score was 43 (SD 16), the median SOFA score was 6.5 (4.0-8.75) and the median NEMS score was 36 (19.25-39.0). All demographic, clinical and laboratory patient characteristics are reported in Table 1.
Regarding the arterial blood gas analyses, the mean pH was 7.39 (SD 0.08), with a mean PaCO 2 of 48 mmHg (SD 12 mmHg), a mean PaO 2 of 82 mmHg (SD 25 mmHg) and a mean bicarbonate level of 28 mmol/L (SD 5 mmol/L). All patients' blood gas analyses results are reported in the Supplementary Material (Table S1); the results of the blood gas analyses for each patient are reported in Table 2.

p50 Shift over Time
For each patient, a p50 calculation was reported, analyzing the mean (SD), median (IQR) and min/max values (Table S2); for all patients, a Kolmogorov-Smirnov test rejected the hypothesis of normal data distribution. Analyzing the whole ICU population, a significant difference between early p50s and late p50s was found (20.63 ± 2.1 mmHg vs. 18.68 ± 3.3 mmHg, t = 2.15, df = 31, p = 0.03), identifying a left p50 shift over time (Figures 2 and 3).  After stratifying patients according to ICU survival, we performed the same analysis; both groups presented a progressive left p50 shift over time. While the early p50s comparison between alive (n = 663 blood gas analyses) and deceased (n = 261 blood gas analyses) patients did not show any difference (t = −2.703, df = 30, p = 0.09), a significant difference was found comparing late p50s between the group of alive patients and the group of deceased patients. Indeed, the p50s of patients who survived were more left shifted than those who deceased (t = −3.29, df = 30, p = 0.01, Figures 4-6). Moreover, matching the median of the p50 in alive compared with deceased patients, a significant difference was found (t = −4.08, df = 7.238, p = 0.004).  . Early/late p50 stratification according to patients' outcome. Boxplot analysis of early/late p50s between patients discharged from ICU alive (n = 663) and deceased (n = 261). To eliminate all confounding factors on p50 modification, such as pH, PaCO 2 , 2,3-BPG and temperature, all p50s were calculated according to the Hill formula, modified by Dash. The circles represent outliers patients. Figure 6. Early/late ODC p50 comparison according to patients outcome. Graphic representation of Hb-O 2 affinity ODC, reporting the curve identified by the median of early p50s (A) and late p50s (B) of patients discharged alive (blue curve) and deceased (red curve). To eliminate all confounding factors on p50 modification, such as pH, PaCO 2 , 2,3-BPG and temperature, all p50s were calculated according to the Hill formula, modified by Dash. We may observe that patients with bad prognosis do not present a significant left shift over time.

Subgroup Analysis
After analyzing the single patients' temporal p50 shift, two different temporal patterns were identified (Figures S3 and S4). The 1st pattern showed a left p50 shift over time and was composed of 22 (68.7%) patients; conversely, the 2nd group presented a relative right p50 shift over time between the early and late phases and was composed of 10 (32.3%) patients. A chi-square analysis did not find any correlation between the pattern of p50  (Table S3).

Variance Analysis
During the ICU stay, p50 progression displayed several variations over time ( Figure  S5). Analyzing the variance of p50s for each patient, important temporal variability was confirmed, ranging from σ 2 = 0.38 to σ 2 = 31.86 (Table S2, Figures S5 and S6). This variability did not result in constant and regular fluctuations during the entire ICU stay but instead showed important irregular fluctuations over time ( Figure S5). A one-way ANOVA performed on each patient's 3rd day groups of blood gas analyses confirmed a p50 temporal variability greater in the early p50s (σ 2 = 8.6, 0.38-31.86) than in the late p50s intervals (σ 2 = 3.84, 0.41-15.12); moreover, a p50 variance reduction over time was found (t = 4.198, df = 31, p < 0.001).

Discussion
Recent literature suggests that critically ill patients admitted for COVID-19 disease present a left shift in ODC, when compared with patients presenting with acute hypoxemia related to other diseases [20]. The present fact may contribute to understand the effects of this new disease on the Hb-O 2 affinity.
In our retrospective analysis, we focused on p50 temporal shift in critically ill COVID-19 patients and found a significant difference in Hb-O 2 affinity at ICU admission compared with the last three days of ICU stay. The identification of a disparity between early p50 and late p50 results suggests that one or more factors modulating Hb-O 2 affinity are involved at these two time points. As the ODC was normalized according to pH, PaCO 2 and temperature, it is relevant to consider that the difference in p50, which appears minimal in our finding, may be relevant from a pathophysiologic point of view. We can in fact affirm that, "despite" the normalization induced by the Hill formula, modified according to Dash, some additional-and thus far, undetermined-factor is able to induce a difference in p50s between ICU admission (early p50) and discharge (late p50).
As already stated by Vogel et al. [20], in our study, the first phase of the disease was characterized by a left shift in ODC, that goes beyond normal values associated with the Bohr/Haldane effect [5]. Therefore, although the ODC in COVID-19 patients is further shifted to the left compared with the standard condition, it showed an interesting tendency toward an absence of left shift in patients with negative outcomes, compared with alive patients. The present trend was not observed in those discharged from the ICU. In other words, a persistent left shift of the curve is observed in patients who will improve their lung function, while, compared with alive patients, an absence of left shift is observed, over time, in patients who will not show a recovery. We can consider, for deceased patients, the present result as a factual absence of p50 left shift trend, taking into account that it was calculated using the Hill coefficient, modified according to Dash (which permits us to avoid any confounding interactions). Based on previous studies, we can assume that, in the early phase, patients' clinical status and lung injury are mainly determined by SARS-CoV-2 infection [29], while, in the late phase, lung inflammation's impact on lungs and on the consequent gas exchange could be assumed to be dissimilar between patients [24,30].
Having assessed this ODC propensity to vary during the clinical course, a COVID-19-related effect (direct or indirect) on ODC appears possible. Some groups, such as Harutyunyan et al. [31] and Wenzhong et al. [18], have suggested that there is a hemoglobin allosteric modulation induced by SARS-CoV-2, possibly due to virally derived, transcribed nonstructural proteins (such as ORF) and surface proteins that can directly bind porphyrins and the heme group [18,32]. As a consequence, we observed a deviation of the ODC out of the normal range, even with the use of the Hill formula modified by Dash, which would translate into a possible change in the hemoglobin quaternary state [31,33]. In clinical practice, the mechanism involved in oxygen transport and release in critically ill COVID-19 patients remains to be deeply investigated [34,35]. However, we may speculate that an overall decreased left shift of the curve over time, in the history of the disease, is probably an indication of aggravation of the respiratory function and consequently of gas exchange. Contrarily, a relative steady left shift of the curve may indicate that the lung status is improving and therefore that tissues are relatively well oxygenated.
This study was burdened by some limits. First, we did not directly measure the viral load in the bloodstream; therefore, we can only hypothesize that it correlates with the p50 change. Studies comparing the viral load in the critical care setting and p50 are therefore needed to confirm this hypothesis; however, this is potentially supported by previous studies [24,30], which have confirmed a higher SARS-CoV-2 viral load at disease onset than at hospital discharge. Second, plasmatic 2,3-BPG measurements were not performed, thus leading to set the value at a default 4.65 mmol/L in all analyzed ABG analyses; although the 2,3-BPG level may have an impact on the p50 shift, in the acute setting, it would not seem to significantly influence the Hb-O 2 affinity [36]. Third, the early and late ABG analyses, grouping in a 3-day interval was arbitrarily defined. The choice was a compromise between having an adequate number of ABG analyses, thus allowing a valid median p50 calculation, and the avoidance of combining too many days together, to be as precise as possible in describing the possible variability and the clinical course. Finally, ICU admission discharge times may not correspond to the same phase of illness for different patients. Although these data must be evaluated cautiously, they allow us to determine both the ODC and p50 shift over time, determining that this shift progressively reduces its fluctuation during the ICU stay.

Conclusions
An important variability and temporal left shift of the p50, ODC and Hb-O 2 affinity in critically ill COVID-19 patients was observed. These data may bring forward an additional step toward understanding the pathophysiology of hypoxemia observed in COVID-19 patients. Interestingly, the absence of a steady left p50 shift over time may correlate with worse clinical outcomes in critically ill COVID 19 patients.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/jcm11030788/s1, Table S1: patient's clinical, laboratory and arterial blood gas analysis characteristics. Table S2: patients' p50 calculated according to Hill formula modified by Dash et al., Table S3: MANOVA analysis comparing p50 shift over time. Figure S1: polynomial expressions from Dash et al., Figure S2: graphic representation of p50. Figure S3: means of early and late p50s for each critically ill COVID-19 patient. Figure S4: of early and late p50s stratified according to ICU outcome. Figure S5: example of temporal evolution of p50 during ICU stay. Figure S6: temporal evolution of p50 variance over days during ICU stay. Figure S7: graphic representation of Hb-O 2 ODC affinity of late p50s in alive patients. Figure S8: graphic representation of Hb-O 2 ODC affinity of late p50s in deceased patients.