The systemic capillaries were modeled as a single compartment (Supplementary Figure S1). All oxygen transfer to the tissues in the systemic circulation was assumed to occur in the capillary region. This was parameterized using an overall metabolic rate for tissue oxygen consumption. Therefore, the total oxygen concentration in the pulmonary arteries was taken to be equivalent to that in the systemic veins ( C p a , O 2 T ≈ C s v , O 2 T ) . Similarly, the total oxygen concentration in the pulmonary veins was taken to be equivalent to that in the systemic arteries ( C p v , O 2 T ≈ C s a , O 2 T ) , assuming that transfer out of the arterioles is not significant. If necessary, oxygen leakage in the arteriole system can be incorporated by equating the total oxygen concentration entering the systemic capillary compartment to a fraction of that leaving the pulmonary capillary compartment ( C s a , O 2 T = f r a c t i o n ⋅ C p v , O 2 T ) .

To determine the total oxygen concentration in the pulmonary arteries, an unsteady-state mass balance was performed over a differential control volume moving through the systemic capillary compartment at the same velocity as the surrounding blood: d C c v , T d t = − M R O 2 V sys,cap (1) Eq. 1 was solved over time to determine the relation between the total oxygen concentration in the pulmonary arteries and veins (Supplementary Appendix S1.1).

An important point to highlight is that the total concentration of oxygen in the blood was taken as the sum of the dissolved and hemoglobin-bound components. The dissolved oxygen concentration is necessary to consider as it drives the gradients for mass transfer, and it is what the tissues are exposed to. In other words, oxygen must dissociate from hemoglobin and dissolve into the blood plasma before it can be transferred to the body tissues. C j , O 2 T t = C j , O 2 d t + 4 ⋅ C H b ⋅ S j , O 2 t (2) In Eq. 2, j = pa, pc, pv, sa, and sv for the pulmonary arteries, pulmonary capillary compartment, pulmonary veins, systemic arteries, and systemic veins, respectively.

The relation between hemoglobin saturation and oxygen partial pressure was estimated using a fit to the standard oxygen-hemoglobin dissociation curve (Severinghaus, 1979): S j , O 2 t = 1 + 23,400 m m H g 3 P j , O 2 t 3 + 150 m m H g 2 P j , O 2 t − 1 (3) P j , O 2 t = C j , O 2 d t β p (3a) In Eqs 3, 3a, j is the same as defined in Eq. 2.

The pulmonary capillaries were modeled as a single compartment, assuming no significant regional heterogeneity in lung ventilation/perfusion (Supplementary Figure S2). An unsteady-state mass balance was performed over a control volume moving within the compartment to obtain a time-dependent, spatial oxygen profile: Δ V d C c v , T d t = k p c C A , O 2 t − C c v , d t (4) k p c = D L , O 2 β p Δ z L c (4a) L c = V p c A eff (4b) C A , O 2 t = β p ⋅ P A , O 2 t (4c)

The total oxygen concentration within the control volume was assumed to be spatially uniform. In addition, the mass transfer coefficient was assumed to be uniform along the length of the pulmonary capillary compartment. Pulmonary blood velocity was used to track the movement of the control volume along the pulmonary capillaries and to allow for the effect of perfusion on the mass transfer mechanism (more details provided in Supplementary Appendix S1.2). In addition, the effect of the O ₂ − Hb reaction was incorporated by using Eqs 2, 2 3, to calculate the dissolved oxygen concentration. The oxygen mass transfer rate was determined using the molar flux across the alveolar-capillary membrane. Pseudo-steady state was assumed to be valid for this transfer process as the diffusion time was much smaller than the time for boundary changes (Supplementary Appendix S3). The dissolved oxygen in the alveolar membrane was assumed to be uniform along the length of the pulmonary capillary compartment. Furthermore, the solubility factor for oxygen was assumed to be the same for the membrane space and blood plasma.

Since the plasma and red blood cells are not separately modeled, the individual diffusion and reaction coefficients (consisting of diffusion through the alveolar-capillary membrane, blood plasma, red blood cells, and facilitated diffusion due to the oxygen-hemoglobin interaction) are not known. Therefore, a physiologic estimate of the lung diffusive capacity was used, which inherently incorporates these terms and is around 21 mL O ₂ ⋅min⁻¹ ⋅mmHg⁻¹ for a normal subject at rest (Guyton and Hall, 2000). Data from a clinical paper showed that the average lung oxygen diffusive capacity did not significantly deviate from this normal value under the condition of anoxia (Fishman, 1954). Therefore, it was assumed to be constant for all cases (simulated and patient).

The alveoli were modeled as a single, homogeneous compartment (Supplementary Figure S1). An unsteady-state mass balance was performed for the inspiration and expiration processes to determine the alveolar oxygen partial pressure in response to input breathing data. The inspired air was assumed to be ideal, saturated with water vapor, at normal body temperature, and at a 0.21 oxygen mole fraction. With no spatial variations in gas concentration, the oxygen partial pressure leaving the alveolar compartment was taken to be equivalent to that within it. Mass transfer out of the alveolar compartment was calculated using the concentration gradient in dissolved oxygen across the alveolar-capillary membrane. A spatially averaged oxygen concentration along the pulmonary capillary compartment ( C c , O 2 ) was used for this gradient. The resulting differential equations are similar to those presented in previous literature (Reynolds et al., 2010).

During inspiration, when d V A d t > 0 , the alveolar partial pressure was determined using: V A t d P A , O 2 d t = d V A d t P I , O 2 − P A , O 2 t − k l C A , O 2 t − C c , O 2 t R ⋅ T (5) P I , O 2 = y O 2 P B − P H 2 O (5a) k l = D L , O 2 β p (5b) C c , O 2 t = 1 L c ∫ 0 L c C p c , O 2 d z , t d z (5c)

During expiration, when d V A d t < 0 , the alveolar partial pressure was determined using: V A t d P A , O 2 d t = − k l C A , O 2 t − C c , O 2 t R ⋅ T (6)

An approximation of a normal breathing pattern was used to validate the model (Reynolds et al., 2010): V A t = 1 2 V vent ⁡ sin 2 π ⋅ b r ⋅ t − π 2 + 1 2 V vent + V End (7) V vent = V T − V D (7a)

Variations of this breathing pattern were used to simulate multiple cases of OSA for assessing the AHI. This was done by varying the breathing rate and tidal volume and introducing apnea/hypopnea events. The AHI for each OSA simulation was defined as the number of apnea/hypopnea events divided by the total breathing time studied (in hours) after achieving model stability. All simulated breathing patterns represented the alveolar volume, and the duration of inspiration and expiration was assumed to be equivalent for each. In these simulations, a sinusoidal function was used to approximate breathing, but the model has the capability to use any form of input breathing data.

Nasal pressure data from two OSA patients was obtained during multi-hour sleep studies, which were conducted with and without continuous positive airway pressure (CPAP) administration. Human protection: the analysis was done using de-identified datasets generated from studies performed in the UCI Sleep Center. The research was done in compliance with the UCI IRB regulations (UCI IRB #267). For each patient, the recorded pressure during a portion of the study without CPAP was converted to a time-dependent lung volume for implementation as a model input. To achieve this, each pressure signal was first normalized to a mean of zero by subtracting a moving average (taken over 70-s intervals for Patient 1 and 50-s intervals for Patient 2) from the raw signal (Supplementary Appendix S4; Supplementary Figures S4, S5). Following this, for each patient, a portion of the signal identified as normal breathing by a clinician was isolated. The breathing rate, determined from the isolated signal, and tidal volume, approximated using the patient’s height-based ideal body weight, were used to simulate a normal breathing pattern and corresponding flow rate for each case (Supplementary Appendix S4; Supplementary Figures S6, S7). Considering that the body-mass index (BMI) falls within the overweight range for Patient 1 (25 ≤ BMI < 30) and the obese range for Patient 2 (BMI ≥30) (Table 3), average functional residual capacities (FRC) measured in overweight and obese subjects suspected of having OSA (Abdeyrim et al., 2015) were used as the end-expiration volumes in the simulated breathing patterns. To convert the recorded nasal pressure to a nasal flow rate, a fitting parameter for each patient was defined by assuming that the average maximum nasal pressure during the identified normal breathing corresponds to the maximum simulated inspiratory flow for ideal tidal breathing. The calculated fitting parameters represent total nasal conductance, effectively combining air density, areas of flow, kinematic heat ratio, and inlet pressure for a laminar flow relationship (Mansour et al., 2002; Mansour et al., 2003). The generated nasal flow signals were then integrated over time to achieve time-dependent lung volumes, which were used to approximate the patient alveolar volumes as model inputs (Supplementary Appendix S4).

All differential equations were solved using the finite difference method (Supplementary Appendix S1). The initial time steps for the simulations and patient cases were chosen to ensure model stability, convergence, and computational efficiency (Supplementary Appendix S5.4). Furthermore, all breathing patterns were introduced after allowing the model to stabilize (t = 360 s). For the results, all presented dissolved oxygen concentrations are in units of μmol/L of blood. In addition, percent decreases in hemoglobin oxygen saturation and dissolved oxygen concentration were calculated by comparing the average normal values to the minimum points during the period of study (Supplementary Appendix S6.2). To determine the percent reduction in the average mass transfer to the tissues, the difference in systemic arterial and venous dissolved oxygen concentrations was averaged over the breathing pattern time and compared to the normal value (Supplementary Appendix S6.2).

For analysis of the clinical performance of the model, each patient was given a proposed hypoxia burden score for a selected event series (Figure 8 for Patient 1 and Figure 9 for Patient 2) and for the entire period of study. To determine the scores, an average and standard deviation of the dissolved oxygen in the systemic arteries across multiple wake sequences were first calculated for both patients (Supplementary Appendix S6.3). A value for the area between the arterial dissolved oxygen curve for the sequence of analysis and the average arterial dissolved oxygen during wakefulness was calculated and normalized by the time of the analyzed sequence to give the hypoxia burden score. Furthermore, an individual burden score was determined for the hypoxia period over the full analysis sequence for each patient. The hypoxia state was defined as a value of the systemic arterial dissolved oxygen lower than one standard deviation below the average arterial dissolved oxygen during wakefulness. The area deviation for this hypoxia sequence was determined as shown in Supplementary Figure S10 and normalized by the total time spent in hypoxia. When assessing the scores, a higher value indicates more severe hypoxia burden. Using the average dissolved oxygen for wakefulness as the threshold to calculate area deviations, the hypoxia burden score for wakefulness would be approximately equal to 0. All equations and further details of the procedure are provided in Supplementary Appendix S6.3.

Furthermore, a statistical analysis was performed on data from Patient 1 to assess the significance of closely occurring obstructive events and to demonstrate the utility of the model for quantitative data processing. The variables studied were the systemic arterial and venous saturations and the difference in dissolved oxygen concentration across the systemic capillaries. Four separate intervals of time were selected, with each one including two obstructive events (apnea or hypopnea) and the subsequent desaturation periods identified in the clinic. An average of each variable was calculated for all four intervals (Supplementary Appendix S6.5). Interval duration was kept constant at 125 s. Overall averages of each variable, taken over the duration of the entire study, were subtracted from the interval averages to create the data set for analysis (Supplementary Appendix S6.5). p-values were computed using the one-sample t-test in MATLAB (“ttest”). The test assesses the hypothesis that the data comes from a distribution with a mean of zero. The changes in the variables due to obstructive events were determined to be significant at p-values lower than 0.05.

The model was first validated by comparing the results for a simulated normal breathing pattern at rest (Figure 1A) to expected physiological values. The average alveolar oxygen partial pressure was around 99 mmHg (Figure 1B). This is within the range of 98–104 mmHg stated in previous literature (Gardner, 1994; Guyton and Hall, 2000). The fraction of hemoglobin oxygen saturation was predicted to be around 0.98 in the systemic arteries and 0.76 in the systemic veins (Figure 1C), which approximates normal values of 0.97 and 0.75, respectively (Guyton and Hall, 2000). Furthermore, on average, the simulated dissolved oxygen concentration was 139 μM in the systemic arteries and 57 μM in the systemic veins (Figure 1D). When converted to oxygen partial pressures ( P s a , O 2 = 99 mmHg and P s v , O 2 = 41 mmHg), the values fall within expected ranges of ≈ 85–100 mmHg and ≈ 27–45 mmHg, respectively (Guyton and Hall, 2000; van Faassen et al., 2009; Ortiz-Prado et al., 2019).

A severe OSA breathing pattern was simulated to further validate the model and to assess the effect of frequent apneas on oxygenation in the systemic arteries and veins (Figure 2). The simulated pattern consisted of four 40-s apneas within a 3.3-min period (AHI = 72) (Figure 2A), which meets the severe OSA criteria of AHI ≥30. These obstructive events were separated by 10-s periods of hyperventilation (V _T = 1.0 L and b _r = 24 breaths/min) (Figure 2A). The resulting breathing pattern is similar to that in previous literature (Netzer et al., 2001; Cheng and Khoo, 2012). Realistically, this pattern would not be repeated multiple times over an hour as a patient would awaken once the hemoglobin saturation reached critically low values. Importantly, the model predicted a progressive decrease in the minimum hemoglobin oxygen saturation and dissolved oxygen concentration with each apnea (Figures 2B, C). For example, in the systemic arteries, the saturation fraction was around 0.93 after the first apnea, dropped to 0.91 after the second one, and further decreased to around 0.90 after the third apnea (Figure 2B). Another significant result is that the alveolar and end-pulmonary capillary oxygen partial pressures were not equivalent at the end of the studied breathing period (Figure 2D). Furthermore, it was noted that there is a brief time-delay from when each apnea starts to when the drop in oxygenation is felt in the systemic arteries and veins and, by extension, in the body tissues (Figures 2B, C).

Two variations of ventilatory response following an apnea were simulated to assess differences in oxygen recovery. In Simulation 1, the hypothetical patient resumed normal breathing (b _r = 12 breaths/min) following each 20-s apnea (Figure 3A). In Simulation 2, there was a 20-s period of hyperventilation (b _r = 24 breaths/min) after each apnea, followed by a resumption of normal breathing (Figure 3B). Considering that these breathing patterns incorporate two 20-s apneas over a 10-min period (AHI = 12), they could represent a case of mild OSA, which has a criteria of 5 ≤ AHI < 15. Based on these simulations, reoxygenation occurred faster when hyperventilation was initiated after an obstructive event. It took approximately 3.5 min, following each drop, to return to the normal saturation and dissolved oxygen concentration with a normal breathing response (Figures 3C, D). However, this time was reduced to around 0.5 min with hyperventilation included (Figures 3C, D). The drops were identical in both scenarios, but the hemoglobin oxygen saturations and dissolved oxygen concentrations increased past the predicted normal values for a short period during the hyperventilatory response.

The decrease for oxygenation in the systemic vessels was compared for some of the simulated cases to assess the relative importance of the AHI and individual apnea duration in determining the severity of decrease in oxygen levels during OSA (Figures 7A–C). For the cases with a constant AHI, the model predicted a positive linear correlation between the apnea duration and the magnitude of percentage decreases (Figures 7D, E). However, there appeared to be no clear correlation between the AHI and percentage decreases (Figures 7A–C).

To demonstrate the clinical utility of the model, recorded heart rate and approximated lung volume data from the converted nasal pressure of two OSA patients were used as time-dependent inputs (Supplementary Figures S11, S13). The model output hemoglobin oxygen saturation in the systemic vessels, which can be compared to the recorded pulse oximeter data (Supplementary Figures S12A, S14A). In some regions for Patients 1 and 2, the recorded S p , O 2 was lower than the arterial saturation predicted by the model, while it was comparatively higher following a respiratory effort related arousal (RERA) event (Figures 8C, 9C), and, in other regions, its fluctuations appeared to be in phase with those of the venous saturation (Figure 9C). In addition, the solution provided the dissolved oxygen concentration in the systemic vessels to better quantify the hypoxic burden on tissues (Supplementary Figures S12B, S14B; Table 4). This gives valuable insight into how the dissolved oxygen in the systemic arteries progressively decreases with continuously occurring obstructive events (Figure 9D). Furthermore, although Patient 2 had a higher AHI score identified in the clinic (calculated as the number of events divided by the total time asleep in hours), the hypoxic burden of Patient 1 was determined to be greater based on the proposed scores (Table 4). This result can be explained by the longer average hypopnea length for Patient 1 from sleep clinic data (Table 4). In addition, the calculated scores allowed for quantification of the overall burden specifically during hypoxic periods in patients, which is valuable insight not reflected in the AHI. For further interpretation of model results, the quantitative analysis performed on data from Patient 1 showed a statistical significance when assessing the differences between the overall and event interval averages for oxygenation (Table 5).