COVID-19 and silent hypoxemia in a minimal closed-loop model of the respiratory rhythm generator

Silent hypoxemia, or “happy hypoxia,” is a puzzling phenomenon in which patients who have contracted COVID-19 exhibit very low oxygen saturation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {SaO}_2$$\end{document}SaO2 < 80%) but do not experience discomfort in breathing. The mechanism by which this blunted response to hypoxia occurs is unknown. We have previously shown that a computational model of the respiratory neural network (Diekman et al. in J Neurophysiol 118(4):2194–2215, 2017) can be used to test hypotheses focused on changes in chemosensory inputs to the central pattern generator (CPG). We hypothesize that altered chemosensory function at the level of the carotid bodies and/or the nucleus tractus solitarii are responsible for the blunted response to hypoxia. Here, we use our model to explore this hypothesis by altering the properties of the gain function representing oxygen sensing inputs to the CPG. We then vary other parameters in the model and show that oxygen carrying capacity is the most salient factor for producing silent hypoxemia. We call for clinicians to measure hematocrit as a clinical index of altered physiology in response to COVID-19 infection.


Background
The global COVID-19 pandemic led to over 1,003,000 deaths in the USA, and over 6,881,000 worldwide, from its onset in late 2019 through March, 2023 (Johns Hopkins University the medulla oblongata (Simonson et al. 2021).ACE2 is the cellular entry point for SARS-CoV-2 (Yuki et al. 2020).Additionally, recent work has shown that there is a shift in the oxyhemoglobin dissociation curve1 in COVID-19 patients (Vogel et al. 2020;Ceruti et al. 2022).Since carotid body chemoreceptors respond to both low O 2 and high CO 2 , a primary problem in these patients may be dysregulation of these sensors and chemosensory reflexes in general.COVID-19 infection has been shown to increase ACE2 expression, leading to changes in sensitivity to both CO 2 and O 2 ; changes in blood gases lead to a concomitant change in activity within the nucleus tractus solitarii (NTS).Recent work has shown that ACE2 is present within the carotid bodies of humans (Porzionato et al. 2021;Villadiego et al. 2021) and there is evidence of altered chemosensation across multiple systems with SARS-CoV-2 infection (Caretta and Mucignat-Caretta 2022).The absence of dyspnea-even though patients exhibit low oxygen saturation-suggests that changes in carotid body inputs to the NTS are a key feature of SARS-CoV-2 infection.Additionally, there may be changes in NTS activity that contribute to the blunted ventilatory response but this has not yet been reported.

Altered chemosensory function and silent hypoxemia
After four years of the COVID-19 pandemic and ongoing endemic infection, a few key pathophysiologies have become apparent.First, ACE2 expression is correlated with the location and severity of infection (Zou et al. 2020).Because ACE2 is, based on current knowledge, the main vector by which SARS-CoV-2 enters the body's cells, changes in ACE2 expression should have an impact on the severity and time course of COVID-19 symptoms.Second, changes in NTS signaling may play a key role in altering the normal, physiological response to changes in oxygenation during COVID-19, and that information may be carried by the glossopharyngeal nerve (innervating the carotid body) or lung afferents via the vagus nerve.Information sensed at the carotid bodies (and lung interoceptors) ultimately reaches the NTS via the vagus and glossopharyngeal nerves.From the NTS, these signals are distributed to local visceral integration circuits within the medulla, including the cardivascular control regions (rostral and caudal in the ventral medulla) and the preBötzinger complex and associated regions of respiratory control within the brainstem.
Based on the clinical observations reported so far, it appears that there is a change in gain in the pathway from carotid body, to NTS, to the breathing rhythm generator and pattern formation network.These observations in patients have provided the motivation for us to focus on assessing the effect of changes in sensitivity/gain in this signaling pathway.This change in gain may be more prevalent in any one of these circuit elements and further work needs to be done to determine the exact mechanism by which sensitivity of the control circuit is impacted.
Given the low partial pressure of oxygen in arterial blood (P a O 2 ) of patients infected with SARS-CoV-2 virus (Sartini et al. 2020;Chen et al. 2020) and the high expression of ACE2 in the carotid bodies, it is likely that altered chemosensory reflexes play a central role in the symptoms and outcomes seen in COVID-19 patients (Porzionato et al. 2021(Porzionato et al. , 2020)).In light of this data, we hypothesize that altered chemosensory function at the level of the carotid bodies and/or the NTS are responsible for this blunted response to hypoxia.

Overview of our approach
We use a previously published computational model of respiratory control (Diekman et al. 2017) to explore this hypothesis by altering the properties of the gain function representing oxygen sensing inputs to the respiratory central pattern generator (CPG).As reviewed in Sect.4, there are several respiratory control models featuring sensitivity to hypercapnia, which under normal circumstances plays the leading role in regulating breathing effort, and few models based on hypoxia-driven chemosensory feedback.Because hypoxia seems to coexist with normal CO 2 levels in silent hypoxemia, we base our investigation on a closed-loop respiratory control model focused on blood oxygen regulation.The respiratory control model studied in Diekman et al. (2017) has seven dynamical variables: voltage of a central pacemaker cell, together with one fast and one slow gating variable; diaphragm muscle activation; lung volume; partial pressure of O 2 in the lung; and partial pressure of O 2 in the bloodstream.Regulation of the endogenous breathing rhythm occurs through hypoxiasensitive chemosensory feedback in the model.Thus, we refer to this system as the 7D-O2 model.The 7D-O2 model strikes a balance between simplicity, in order to preserve analytic transparency, and complexity, in order to capture the phenomenon of interest.See Sect. 4 for further discussion about the realism/tractability trade-off in modeling respiratory control.
In this paper, we use the 7D-O2 model to explore our hypothesis by altering the properties of the gain function representing oxygen sensing inputs to the CPG.We then vary other parameters in the model, and show that oxygen carrying capacity is the most salient factor for producing silent hypoxemia.We exploit the presence of a small parameter (the Henry's Law constant) in the expression for the O 2 saturation curve to provide a mathematical explanation for the effect of Fig. 1 Schematic of the 7D-O2 model.Bursts of action potential firing (V , mV) in the respiratory central pattern generator (CPG) drive a pool of motor neurons (α, dimensionless), leading to expansions of lung volume (vol L , L) and increases in lung and blood oxygen (P A O 2 and P a O 2 , mmHg).Through a chemosensory pathway (g tonic , nS), the blood oxygen level affects the amount of excitatory current sent to the CPG, thereby closing the control loop (red arrow).Time (t, seconds) is the horizontal axis for all traces.The seven parameters shown in blue (chemosensory feedback parameters φ, θ g , and σ g , see Methods; hemoglobin concentration [Hb]; base lung volume vol 0 ; time constant τ LB ; and metabolic demand parameter M) are varied in this study to model silent hypoxemia.Redrawn, with modifications, from Diekman et al. (2017) changing the hemoglobin concentration (hematocrit) in the model. 2.
Figure 1 shows a schematic of the 7D-O2 model with components representing CPG membrane potential (V ), motor pool activity (α), lung volume (vol L ), lung oxygen (P A O 2 ), blood oxygen (P a O 2 ), and chemosensation (g tonic ).The model has a closed-loop structure since an excitatory current, I tonic , depends on P a O 2 and is an input to the CPG component (red arrow).In this model, the rate of metabolic demand for oxygen from the tissues is represented by the parameter M. If metabolic demand is low or moderate (M < 1.2×10 −5 ms −1 ), then the model exhibits a stable eupneic rhythm with CPG bursting activity driving fluctuations in lung volume that bring in a sufficient amount of oxygen to maintain P a O 2 in the normoxia range (see the "plateau" region of the P a O 2 versus M curve shown in Fig. 2a and the traces in the left panel of Fig. 2b.)However, if metabolic demand is too high (M > 1.2 × 10 −5 ), then the model exhibits a form of tachypnea that transitions to tonic CPG bursting activity that does not drive the lungs to effectively maintain P a O 2 in the normoxia range (see the "collapse" region of Fig. 2a and the right panel of Fig. 2b).
In silent hypoxemia, we would expect to observe a lower height for the plateau region of the P a O 2 versus M curve, since these patients display abnormally low P a O 2 despite minimal changes in minute ventilation.There are three possibilities regarding the collapse region in silent hypoxemia patients: the collapse point could shift to a lower M value (as illustrated in Fig. 2c), stay at the same M value (as in Fig. 2d), or shift to a higher M value (as in Fig. 2e).Due to an acute disease-induced reduction in steady-state P a O 2 , it seems plausible that there would be a decrease in tolerance of higher metabolic demand.Therefore, we explored parameter space to see if the closed-loop model is capable of producing P a O 2 versus M curves with shapes similar to the hypothetical curve shown in Fig. 2c.

Quantitative modeling approach
Quantitative modeling has helped elucidate principles of normal and pathological functioning of the respiratory system, although its fundamental mechanisms remain debated.Mathematical models can be particularly helpful for generating experimentally testable hypotheses.A variety of models have been developed for the respiratory CPG (Butera et al. 1999a, b;Del Negro et al. 2002;Rubin 2008;Del Negro and Hayes 2008;Rubin et al. 2009;Phillips and Rubin 2019;Phillips et al. 2022), for chemosensory feedback-based regulation schemes (Grodins et al. 1954;Khoo et al. 1982;Cherniack and Longobardo 2006), and for cardiopulmonary gas exchange (Ben-Tal 2006).See Molkov et al. (2017) and Lindsey et al. (2012) for a review.A smaller number of published models represent closed-loop control incorporating a conductance-based CPG, muscle dynamics, gas exchange, and sensory feedback (Ben-Tal and Smith 2008;Park et al. 2012;Molkov et al. 2014).Of these, several focus on hypercapnia (excessive CO 2 ) as the regulatory pathway.In order to generate hypotheses about silent hypoxemia, we chose to work with a conductance-based CPG model with O 2 chemosensation as the sensory feedback pathway closing the control loop.To our knowledge, our previously published model (Diekman et al. 2017) is the only model meeting these criteria.Aspects of it have been experimentally validated (Diekman et al 2022;Diekman et al. 2018).Like any model, this model fails to represent all aspects of the control system.We have not included CO 2 sensing in our model due to the high diffusion rates of CO 2 when compared to O 2 in the lung (West 2008) and evidence showing that CO 2 is ≤ 35 mmHg in patients presenting with silent hypoxemia (SH) and minimal tachypnea (Chandra et al. 2020;Alamé et al. 2022).Additionally, we do not explicitly include rapidly adapting (RAR) or slowly adapting (SAR) lung mechanoreceptors in the model-lung volume is present in the model and repro-duces inspiratory drive in much the same way that SARs do in vivo.Nevertheless, in spite of these limitations, the model suffices to generate testable hypotheses that could be pursued by the clinical community.
The 7D-O2 model is a closed-loop respiratory control model that comprises a well-established conductance-based central rhythm generator (the Butera-Rinzel-Smith model (Butera et al. 1999a;Diekman et al. 2017)) with a voltage variable V , a fast gating variable (delayed-rectifier potassium current activation, n), and a slow gating variable (persistent sodium current inactivation, h).The output of the BRS model cell, namely the voltage, drives a motor pool activation variable, α, that in turn drives expansion of the lungs.The lung volume (vol L ), the partial pressure of oxygen in the lungs (alveolar pressure, P A O 2 ), and the partial pressure of oxygen in the bloodstream (P a O 2 ) complete the model variables.The BRS cell includes an excitatory current driven by a tonic conductance that is regulated by chemosensory feedback, closing the control loop.When the tonic conductance assumes intermediate values, the BRS cell exhibits bursting activity, consistent with eupnea (normal steady breathing).If blood O 2 levels are significantly reduced, the tonic conductance increases, which can trigger a transition into a rapid, tonically firing "beating" regime, analogous to tachypnea (pathologically rapid shallow breathing).If blood O 2 levels are significantly increased, the tonic conductance decreases, which can push the BRS model cell into a stable resting fixed point at, which Butera et al. called the "quiescent" regime (Butera et al. 1999a).The 7D-O2 model includes a metabolic demand parameter, M, regulating the rate at which oxygen is removed from the bloodstream to the tissues.As the "phenotype" or "physiology" of the model, we take the steady-state value of P a O 2 as a function of M. For the original model as presented in Diekman et al. (2017), the P a O 2 -vs-M curve shows a plateau near 100 mm Hg (normoxia) that collapses to a critically hypoxic state when M increases past a high threshold (Fig. 2a).As we varied the original parameters to investigate possible mechanisms of silent hypoxemia, we monitored the height of the normoxia plateau, and the location of the collapse point.
Simulations were conducted using MATLAB v. 2020B on the NJIT high-performance computing cluster.Code corresponding to the 7D-O2 model is posted on Github at https:// github.com/ModelDBRepository/229640and on ModelDB at https://modeldb.science/229640.Code corresponding to the silent hypoxemia model developed here is available at https://modeldb.science/2015954.

Relating model parameters to potential silent hypoxemia mechanisms
The mechanism by which COVID-19 leads to sustained hypoxemia in the absence of dyspnea is currently unknown.
The minimalist model of Diekman et al. (2017) includes a number of key parameters that are plausible targets for modification to mimic the effects of COVID-19-infection on respiratory control.Oxygen carrying capacity is a key variable in pulmonary mechanics.Repeated bouts of intermittent hypoxia, as seen in obstructive sleep apnea, can increase HIF-1α signaling, with a subsequent increase in erythropoietin, and an increase in hemoglobin and erythrocytes.Similar changes are seen in conditions that result in chronic hypoxemia and hypercapnia, such as cardiovascular disease, obstructive sleep apnea, and chronic obstructive pulmonary disease (Mauad et al. 2021;Paquette et al. 2021;Li et al. 2019).Many of the patients presenting with silent hypoxemia have pre-existing conditions and comorbidities that are likely to increase hematocrit and this increase in oxygen carrying capacity may blunt chemoreceptor responses-exacerbating the "happy hypoxia" phenomenon.Unfortunately, no current literature quantifies hematocrit in these patients.
Motivated by these observations, we systematically varied (plus or minus 20%) the following parameters that control the saturating effect of hypoxia-sensitive chemosensory feedback to the central pattern generator: σ g , which controls the slope of the sensory feedback curve at maximum sensitivity (gain at threshold); θ g , which controls the threshold activation value for sensory feedback (50% activation point); and φ, which controls the maximum sensory feedback drive at full activation.
Lung volume is a key determinant of mechanosensory feedback to the NTS and the CPG.Our model incorporates lung volume and allows us to monitor changes in lung volume in response to changes in central drive for breathing.This also allows us to monitor lung volume as an outcome measure to determine if the CPG is actually causing lung inflation in a way that assures sufficient gas exchange to sustain life when extrapolated to animal models or human subjects.
Ventilation/perfusion matching is a key drive for respiration.In mammals, the interplay between cardiovascular and respiratory control is essential for ensuring that sufficient oxygen is delivered to the body and CO 2 is removed via the lung.We have included a time constant for O 2 transport between the lung and blood which allows us to simulate changes in diffusion and dwell time within the lung that correlate with diseases such as chronic obstructive pulmonary disease (COPD) and lung fibrosis.Oxygen consumption and CO 2 production are key elements for determining how changes in breathing can match metabolic demand.We have included a simplified treatment of metabolism in the model.As a "biomarker" to test the model behavior, for all parameter sets we varied the metabolic demand parameter M across a range of values.We have not included CO 2 in this model, because CO 2 diffuses up to 20 times faster than O 2 (West 2008) and patients with SH do not appear to be hypercap-nic since there is very little change in breathing rate-CO 2 is a potent stimulator of minute ventilation and hypercapnia results in pronounced increases in breathing frequency (Moosavi et al. 2003;Parshall et al. 2012;Nakano et al. 2015).
Additionally, we vary the hemoglobin concentration to mimic the effect of chronic hypoxia seen in humans living in hypoxic environments which can include mountain dwellers (Hancco et al. 2020), individuals with severe obstructive sleep apnea (Li et al. 2019), or other cardio-respiratory disorders (Paquette et al. 2021;Balasubramanian et al. 2021).These individuals can have high hematocrit, a corresponding increase in red blood cells, and increased blood viscositysimilar to what has been reported in COVID-19 patients (Choi et al. 2022).

Twenty percent variation in parameters specifying chemosensory feedback gain is insufficient to qualitatively reproduce silent hypoxemia
Motivated by the hypothesis that silent hypoxemia could result from a dysregulation of carotid body O 2 receptors, we first considered variation of the parameters associated with the chemosensory pathway of the model.In the 7D-O2 model, there is a sigmoidal relationship between P a O 2 and g tonic , with the parameters φ, θ g , and σ g controlling the height, half-activation, and slope of the sigmoid, respectively (see Fig. 3a).We simulated the closed-loop model over a range of M values while varying these parameters over three levels spanning roughly ±20% of their original values (φ = 0.24, 0.3, 0.36, θ g = 70, 85, 100, and σ g = 0.24, 0.3, 0.36), yielding 27 different combinations in total.Figure 3b shows that varying these parameters generates P a O 2 vs M curves in which the plateau and collapse point are shifted down and to the right (similar to the hypothetical case shown in Fig. 2e).None of the 27 combinations, however, produce any curves with the plateau and the collapse point shifted down and to the left (similar to Fig. 2c).In patients with comorbidities that cause compensatory adaptations to chronic hypoxiadownstream from hypoxia-inducible factor 1α (HIF-1α)-it seems plausible that disease-induced reduction in steadystate P a O 2 would be accompanied by an increase in tolerance of higher metabolic demand.However, we do not consider any of these model variants to suitably capture the phenomenon of silent hypoxemia.
In order to proceed further, we selected a single parameter set from among the 27 combinations previously considered as our working model for producing the hypoxic plateau region, namely φ = 0.24, θ g = 70, and σ g = 36.These parameters gave the curve with the greatest reduction of P a O 2 (dark-est blue curves in Fig. 3a, b), although the collapse point did shift to significantly higher values of M. Next, based on reports indicating that COVID-19 patients have altered oxyhemoglobin dissociation curves (Vogel et al. 2020;Ceruti et al. 2022), we considered variation of the model parameter K which represents hemoglobin binding affinity (Eq.18 in Methods).

Varying the shape of the hemoglobin saturation curve leaves blood oxygen unchanged and weakly shifts the metabolic collapse point
The effect that increasing the binding affinity (decreasing K ) has on the hemoglobin saturation versus blood oxygen curve (SaO 2 − P a O 2 ) with the new chemosensory parameters is shown in Fig. 4a.Tighter binding affinities (K values less than the default value of 26 mmHg) do shift the P a O 2 vs M curve to the left, but the respiratory collapse point is still at higher metabolic demand values than the original model (Fig. 4b).See also Fig. 9a in "Appendix A.2" for the effect of varying K with the original chemosensory parameters.
Since varying the chemosensory parameters alone was not sufficient to model a silent hypoxemia patient prone to respiratory collapse, we considered other parameters that could plausibly be affected by COVID-19.Specifically, lung damage due to excessive immune response or local thrombosis could reduce the effective unloaded lung volume (model parameter vol 0 ), or impede the flux of oxygen between the alveoli and the alveolar capillaries.The latter effect could be reflected by an increase in the model parameter τ LB , which governs the effective relaxation time for differences in partial pressure of oxygen in the model's lung and blood compartments, respectively.Therefore, while keeping the chemosensory sigmoid parameters (φ = 0.24, θ g = 70, σ g = 36) fixed, we varied the unloaded lung volume (vol 0 = 1.6, 2.0, 2.4) and the time constant for the flux of oxygen from the lung to the blood (τ LB = 100, 500, 900).demand.The problem remains of finding parameters that can shift the collapse point without elevating the eupneic plateau.

Varying both oxygen flux and lung volume has little effect on the blood oxygen vs metabolic demand curve
Another parameter that could possibly be affected by COVID-19 infection is hematocrit (hemoglobin concentration).Increased hematocrit (polycythemia) is one phenotypic response observed in individuals who relocate from sea level to extreme high altitude environments for a prolonged period of time (Winslow and Cassinelli 1987;Beall et al. 1990).

Increasing hemoglobin concentration shifts the P a O 2 collapse point to lower M values while maintaining eupneic plateau height
Finally, we considered variation of the parameter [Hb] representing the hematocrit, i.e., the concentration of hemoglobin within the blood, which was set to 150 g/l in the original 7D-O2 model.Figure 5b shows that increasing [Hb] within the model lowers the collapse threshold of the P a O 2 versus M curve, while maintaining a hypoxemic plateau around 80 mmHg.A 33% increase in [Hb] shifts the collapse point to a similar M value as the original 7D-O2 model, consistent with the hypothetical silent hypoxemia P a O 2 vs M curve shown in Thus, we will consider the model with [Hb]=250 (the second curve from the left in Fig. 5b) as a our working model for silent hypoxemia, and analyzed the model dynamics for simulations in the plateau region and in response to increases in metabolic demand.
Figure 6 compares simulations of the silent hypoxemia model (blue traces) and the original 7D-O2 normoxia model (red traces) for different values of the metabolic demand parameter M, generated as follows.First, the corresponding model is simulated with M = 0.8 × 10 −5 ms −1 for two minutes of simulated time, to establish baseline initial conditions on the eupneic limit cycle.Then, the value of M is changed to the value shown on the horizontal axis (central column: panels b, e, h, k) and the simulation is run for another 10 min of simulated time.For the conditions shown in detail in the left column (M = 0.4 × 10 −5 ms −1 , panels a, d, g, j), this duration is sufficient to effectively remove tran-sient behavior.For the conditions shown in the right column (M = 0.97×10 −5 ms −1 , panels c, f, i, l), the transient effects are still visible for the silent hypoxemia curves. 3 Figure 6a shows voltage traces in the plateau region (M = 0.4×10 −5 ms −1 ) for both the silent hypoxemia model (blue) and the original 7D-O2 normoxia model (red).The frequency of bursting is similar in the two models (Fig. 6b), but there are a few more spikes per burst in the hypoxemia model (Fig. 6d,e).This leads to slightly more vigorous lung expansions in the hypoxemia model (Fig. 6g); however, the levels of oxygen in the blood remain substantially lower (Fig. 6j).As the metabolic demand is increased, the frequency of bursting in the hypoxemia model becomes much faster than in the normoxia model (Fig. 6b), and there are substantially fewer spikes per burst (Fig. 6e).This type of bursting activity leads to more frequent but less vigorous lung expansions and ultimately respiratory collapse at lower levels of metabolic demand in the hypoxemia model compared to the normoxia model (Fig. 6e,k).For comparison, panels (c, f) show voltage traces for a metabolic demand value M = 0.97 × 10 −5 ms −1 for which the SH model's breathing pattern has moved past the point of collapse, while the original 7D-O2 model maintains steady eupneic breathing; note the dramatically reduced mean P a O 2 values in panel (l).
Figure 6h plots minute ventilation (MV) as a function of metabolic demand M for the original 7D-O2 normoxia model (red) and the silent hypoxemia model (blue).MV is a well-established clinical measure of respiratory performance, and is defined as the net volume of respired air.MV is approximately six liters per minute in normal, resting adults, and typically increases with modestly increasing metabolic demand.For excessively high demand, MV shows nonmonotonic behavior in our model, first increasing and then rapidly decreasing.For our model systems, we define MV as the net inspired air per breath (maximum lung volume minus minimum lung volume), divided by the breath cycle duration.For the original model parameters, MV begins near 3 l/min at low metabolic effort (M = 0.2×10 −5 ms −1 ), and increases gradually to approximately 12 l/min at intermediate effort (M 1.2×10 −5 ms −1 ) before collapsing to near zero at excessively high effort (M 1.2 × 10 −5 ms −1 ).In contrast, for the silent hypoxemia model, MV begins with slightly elevated values, relative to the normoxia model, climbs gradually while remaining slightly above the normoxia curve, until suddenly collapsing at MV ≈ 1.0 × 10 −5 ms −1 .For comparison, panels (I,L) shows lung volume and P a O 2 , respectively, for the higher demand value M = 0.97 × 10 −5 ms −1 .The traces show the ongoing transient decline of the respiratory pattern in the SH model, toward full tachypneic collapse, after ten minutes of elevated metabolic demand. 3The same procedure was used to generate Figs.3b, 4b, and 5.
Having established our working model for silent hypoxemia, we next exploit the relative simplicity of the model to investigate the underlying mechanism by which changing hematocrit shifts the collapse point along the M-axis.

Dimension reduction via fast-slow analysis shows varying hematocrit levels has similar effects in the model
Fast-slow dissection is a principled approach to understanding the behavior of dynamical systems involving variables with disparate timescales (Fenichel 1979;Rubin and Terman 2002).Let the vector x represent the fast variables and let the scalar y be the slow variable in the two-timescale system where 0 is a small parameter.Suppose (1), the fast subsystem, has either a stable fixed point x = x fp (y) or else a stable limit cycle solution with period T (y), which we write as x = γ (y, t) = γ (y, t + T (y)), for each value of y in the relevant range.The value of the fixed point, or the shape and period of the limit cycle trajectory, may depend on y.In the limit of small , the variable y given by (2), the slow subsystem, is approximately constant.Then (2) may be written in terms of rescaled (slow) time τ = t as where ḡ(y) represents the average effect of the fast subsystem on the slow subsystem, at a given value of the slow variable y: When the timescales of y and x are sufficiently separated, the dynamics given by (3) provide a lower dimensional approximation of the full system (1)-(2).We previously determined P a O 2 to be the slowest dynamical variable in the 7D-O2 model (Diekman et al. 2017).Identifying the slow variable y with P a O 2 , we applied fast-slow dissection as described above.
Figure 7 shows the averaged rate of change of the slow subsystem, ḡ, defined by Eq. ( 4  The curves in panels a, b have a detailed structure related to the bursting, beating and quiescent regimes of the original BRS model (Butera et al. 1999a).When P a O 2 is in an intermediate range, roughly 70-90 mm Hg in the approximate reduced model, the fast subsystem (1) is in the bursting regime, which efficiently drives gas exchange in the lungs, so that more O2 enters the blood stream than leaves it ( ḡ > 0) provided M is not too large (green and blue curves, Fig. 7a).Within the eupneic range, the ḡ versus P a O 2 curve has a scalloped shape due to the addition of spikes to the bursting pattern as P a O 2 increases.When P a O 2 is below the eupneic range, g tonic increases, forcing the fast subsystem into the steady spiking or "beating" regime, which leads to inefficient gas exchange in the lungs and low minute ventilation.Under these conditions, less O2 enters the blood stream than leaves it ( ḡ < 0).When P a O 2 is above the eupneic range, i.e., P a O 2 90 mm Hg, g tonic decreases sufficiently that the fast subsystem enters the quiescent state, as described in Butera et al. (1999a).That is, the voltage and gating vari-ables enter a stable fixed point corresponding to a steady resting potential.Under these circumstances, no new oxygen enters the blood stream; meanwhile, O2 leaves in proportion to P a O 2 , so the level curves of ḡ decrease rapidly.Referring to Eqs. ( 20)-( 25), we see that if the slow variable y = P a O 2 is held constant in a range where the fast subsystem enters the quiescent state then the expression for ḡ simplifies to giving the smooth descending curves above P a O 2 90 mm Hg (see Fig. 7a, arrows).In the second line, we have used the fact that β O 2 1.We discuss this small parameter further below.In the SH model (Panel b), a similar structure is apparent, but is shifted to the left along the P a O 2 axis.
The curves showing when P a O 2 will increase ( ḡ > 0) and decrease ( ḡ < 0) provide a simplified explanation of the collapse from eupnea to tachypnea through a saddlenode bifurcation in the slow subsystem.In Fig. 7a, b, zero-crossings of ḡ with positive and negative slopes correspond to unstable and stable fixed points of the slow subsystem, respectively.For the 7D-O2 model (Panel a), we can observe the following (cf.Diekman et al. (2017)).When M = 0.4 × 10 −5 ms −1 (green curve), the system has a stable fixed point corresponding to eupneic bursting (P a O 2 =89 mmHg), a stable fixed point corresponding to tachypneic spiking (P a O 2 =41 mmHg), and an unstable fixed point (P a O 2 =74 mmHg).When M = 0.8 × 10 −5 ms −1 (blue curve), the system still has two stable fixed points, but the stable eupneic point (P a O 2 =87 mmHg) and the unstable fixed point (P a O 2 =80 mmHg) have moved closer together.Further increases in M lead to a saddle-node bifurcation in which the stable eupneic point and the unstable fixed point collide and disappear, leaving only the tachypneic fixed point.For example, when M = 1.6 × 10 −5 ms −1 (magenta curve), the system has only 1 fixed point, which corresponds to stable tachypneic spiking (P a O 2 =17 mmHg).For the SH model (panel b), as in panel a the system again has three fixed points for M = 0.4 × 10 −5 ms −1 and only 1 fixed point for M = 1.6 × 10 −5 ms −1 ; the qualitative behavior is the same, although the value of M at which the saddle-node bifurcation occurs is different.
Panels c, d of Fig. 7 show ḡ as a function of both blood oxygen (P a O 2 ) and metabolic demand (M), for the reduced system (4) in the 7D-O2 model (panel c) and the SH model (panel d).The black curve in each panel shows the location (P a O 2 value) of fixed points ( ḡ = 0) in the averaged slow subsystem as a function of metabolic demand M. The heatmap colors indicate the value of ḡ.For the 7D-O2 model (panel c), at M = 0.25 × 10 −5 ms −1 , the lower stable branch and unstable middle branch collide and these fixed points are destroyed in a saddle-node bifurcation (SN 1 ), leaving only the stable upper branch for M < SN 1 .Similarly, at M = 0.88 × 10 −5 ms −1 , the upper stable branch and unstable middle branch collide in another saddle-node bifurcation (SN 2 ), leaving only the stable lower branch (tachypnea) for M > SN 2 .Panel d shows qualitatively similar behavior for the SH model, but the "eupneic" region is shifted to lower P a O 2 values, and the ḡ = 0 curve is shifted to higher values of M for corresponding values of [Hb].Panels e, f illustrate how the curve of fixed points shifts to the right as [Hb] is decreased (blue) and to the left as [Hb] is increased (red) for the normoxia (panel e) and silent hypoxemia (panel f) chemosensory parameters.

Why does changing [Hb] shift ḡ along the P a O 2 axis?
Reexamining the detailed model equations ( §A) we see that the metabolic demand parameter M occurs only in equation (23).We further note that equation ( 23) contains a small parameter, namely the Henry's Law constant (β O 2 ) representing the solubility of oxygen in the blood, in the absence of hemoglobin.For the physiologically realistic values chosen in the original model, 3 ≈ β O 2 P a O 2 η SaO 2 ≈ 150, in appropriate units.Given the form of ( 23), it is clear that setting β O 2 ≈ 0 is a small (regular) perturbation of the dynamics.Neglecting this small parameter, we see that the flux of oxygen from the blood to the tissue is mainly driven by the product M × η of the metabolic demand parameter M and the hematocrit (concentration of hemoglobin) parameter η.Therefore, the model dynamics are (approximately) invariant to any rescaling M → γ M, η → η/γ , i.e., any rescaling of M and η that leaves the product Mη constant.We thus may predict that running the model with η increased and M reduced in proportion would shift a P a O 2 -versusln(M) curve to the left with little deformation.Thus, if (in the model) the results of COVID-19 infection lead to an increase in a patient's hematocrit level, for instance, through hypoxiadriven polycythemia, together with the parameter changes necessary to lower the eupneic P a O 2 plateau as in Fig. 5b, we might expect to recover a net shift of the P a O 2 -v-M curve down and to the left, as suggested in Fig. 1 panel e.
It is natural to conjecture that a similar rescaling might be observed in the full system, since the small-parameter argument above applies equally well to the full seven-dimensional model.Figure 8

Discussion
In this study, we applied a previously published model to the problem of silent hypoxemia seen in a subset of COVID-19 patients.While our model is highly simplified, it can nevertheless be useful for testing hypotheses (about the model) and for generating novel hypothesis relevant to the clinic-although (we hasten to note) we do not make recommendations for clinical practice.We consider a simplified breathing model even though more complicated models are available (Molkov et al. 2017).In general, the complexity of a biophysical model should be related to the research question the model is used to address (Thorburn 1918;Koch and Segev 1998;Carnevale and Hines 2006;Levenstein et al. 2023).Simple models often provide useful tools for understanding the behavior of control systems.In the literature, one may find reduced models of the central controller for breathing that are less representative of preBötzinger complex neurons than our version of the Butera-Rinzel-Smith model (Khoo et al. 1982;Khoo 2000;Cheng et al. 2010).Thus, our intention here has been to develop a model that is as reduced as possible, yet still able to recapitulate changes in control of breathing that can be seen in a given disease condition.Our goal has not been to capture the full complexity of the rostroventrolateral medulla and the entirety of the brainstem breathing control network.However, we have included aspects of the control circuitry that are important in any whole-body model of breathing control.These include the closed-loop feedback, which represents chemosensation (of oxygen levels via P a O 2 ) and mechanosensory inputs, in an abstracted form at least, via the lung volume component of the model.Because we see responses to changes in breathing patterns that recapitulate many of the features observed in COVID-19 patients in a clinical setting, we see our model as a useful tool for exploring the overall control system for breathing.Our ultimate goal will be to increase the complexity of our model to more closely resemble the architecture and diversity of the brainstem control circuitry.However, that is beyond the scope of this manuscript.
In order to generate hypotheses about silent hypoxemia, we chose to work with a conductance-based CPG model with O 2 chemosensation as the sensory feedback pathway closing the control loop.To our knowledge, our previously published model (Diekman et al. 2017) is the only model meeting these criteria and our goal here was to extend that model to address a relevant clinical problem.Aspects of our model have been experimentally validated (Diekman et al 2022;Diekman et al. 2018).However, as with any computational model, our model is not designed to encompass all aspects of the respiratory control system.Despite its limitations, the model suffices to generate hypotheses that can be tested in animal models and by the clinical community.
For the original model as presented in Diekman et al. (2017), the P a O 2 -vs-M curve shows a plateau near 100 mm Hg (normoxia) that collapses to a critically hypoxic state when M increases past a high threshold (Fig. 2a).The work we report here focuses on expanding the original parameters to investigate possible mechanisms of SH.Changing these parameters allowed us to monitor the height of the normoxia plateau, and the location of the collapse point.We hypothesized that altered chemosensory input to the carotid bodies and, eventually, to the NTS and the rest of the breathing control circuitry, is a key factor in silent hypoxemia.However, our simulation results suggest that while changes in chemosensivity may play a role in silent hypoxemia, changes in metabolism and oxygen carrying capacity may have greater relevance for replicating the respiratory collapse seen in these patients.Specifically, altered chemosensitivity can create a hypoxemic plateau region (SaO 2 < 90 mmHg) for a broad range of metabolic demand levels (M = 0.4 to 1.5×10 −5 ms −1 , see blue curves in Fig. 3b).When hemoglobin concentration is then increased, moderate levels of metabolic demand (M = 0.8 to 1.0×10 −5 ms −1 ) lead to complete respiratory collapse (SaO 2 < 60 mmHg, see blue and purple curves in Fig. 3f).
Our hypothesis was based on the premise that O 2 sensing is the key factor in SH.Canonically, it has been suggested that CO 2 is a primary driver for dyspnea (Cherniack and Altose 1987;Chonan et al. 1990;Guyenet and Bayliss 2022), but there is evidence that both hypoxia and hypercapnia equivalently drive the sensation of air hunger (Moosavi et al. 2003).However, clinical case and cohort studies show that patients with SH are not hypercapnic (Chandra et al. 2020;Alamé et al. 2022).This suggested to us that dysregulation of O 2 sensation is a key contributor to the issues seen in SH.We tested this hypothesis by changing O 2 sensitivity in the model at the level of the carotid bodies and NTS as well as evaluating whether those changes could reproduce the SH phenotype.
Complicating factors for these patients include comorbidities that show correlation with poor outcome in patients with COVID-19.These comorbidities include obstructive sleep apnea (OSA), chronic obstructive pulmonary disease (COPD), cardiovascular disease (including hypertension or heart failure) among others.Patients suffering from these diseases often develop polycythemia-an increase in the hemoglobin and hematocrit to adaptively increase the O 2 carrying capacity of the blood.High-altitude populations are well adapted to chronic hypoxia and typically have a higher hematocrit in Andean populations versus Himalayan high-altitude dwellers (Beall and Reichsman 1984), likely due to different adaptation mechanisms.However, subjects with cardiovascular disease (Valeanu et al. 2022), obstructive apnea (Rha et al. 2022), and familial hyperlipidemia (Paquette et al. 2021) also show increased hematocrit.
One consequence of pumping thicker blood is to increase the metabolic demand, even during rest.As we show in our results, increasing metabolic demand increases the likelihood of respiratory collapse.Somewhat paradoxically, with an increased oxygen carrying capacity, the patient may be less able to compensate for the worsening P a O 2 and a critical tipping point for metabolic demand is reached where respiratory efforts are insufficient to keep up with demand.We have not yet seen any report documenting changes in hematocrit in COVID-19 patients who exhibit silent hypoxemia.Based on our modeling results, we would predict that these patients may show increased hematocrit levels.In support of our prediction, a recently published study (Choi et al. 2022) showed that higher blood viscosity was associated with an increase in mortality in COVID-19 patients.Obtaining this kind of data should be possible for patients admitted to the intensive care unit and should be a priority for future investigation.
Angiotensin-converting enzyme 2 (ACE2) is expressed in the lungs, carotid bodies, and respiratory region of the brainstem, and is likely the vector by which the SARS-CoV-2 virus invades the carotid bodies and/or the NTS, thereby potentially contributing to silent hypoxemia.High ACE2 levels also occur in the most vulnerable target organ systems seen in COVID-19 (elevated expression levels occur in lung, heart, ileum, kidney, and bladder (Zou et al. 2020)).Since ACE2 expression is very high in the lungs, and since diffuse alveolar damage, bronchopneumonia, and alveolar hemorrhage are common in COVID-19 (Mauad et al. 2021), it seems reasonable to hypothesize that the decrease in gas exchange across the alveolar membranes within the lung can alter not just the O 2 carrying capacity, but also increase metabolic demand for perfusion of the damaged lung.It may be of value to assess differences in mitochondrial activity in lung cells from normal and COVID-19 patients, or in animal models that have used SARS-CoV-2 or spike protein (now commercially available) to mimic the lung damage seen in human patients.Such experiments would provide data concerning cellular metabolism and give us greater understanding of the impact COVID-19 has on metabolic demand at all tissue levels.
Lack of dyspnea (breathing discomfort) in patients arriving at already overcrowded emergency rooms leads to triaging patients who are not in obvious respiratory distress, when in fact these patients often have reduced oxygen saturation (Bertran et al. 2020).Perhaps the greatest mystery that remains unresolved is why dyspnea is not typically seen in patients exhibiting silent hypoxemia.Sensory perception is subjective and can vary with a host of factors that include sex, socioeconomic background, and ethnicity (Green et al. 2003;Reynolds Losin et al. 2020;de Araújo Palmeira et al. 2011).There is some controversy about these correlates but they may be underlying factors that influence the reporting of silent hypoxemia.Once again, some demographic data is available concerning COVID-19 infection, mortality, and morbidity, but this information has not yet been correlated with silent hypoxemia.Ideally, demographic factors should be reported along with other patient data to better understand the incidence and severity of silent hypoxemia and dyspnea.
Patients with COVID-19 are also subject to mitochondrial dysregulation that contributes to severity and lethality.Mitochondrial function is impacted by the "cytokine storm," a hallmark of the immune response to COVID-19.Thus, upregulation of cytokine release in the context of comorbidities that increase inflammation, including metabolic syndrome, obesity, type 2 diabetes, and increasing age-in addition to the lung and cardiovascular diseases mentioned above, are all associated with mitochondrial dysfunction (Moreno Fernández-Ayala et al. 2020;Grossini et al. 2021;Wang et al. 2022).SARS-CoV-2 infection causes multi-system changes at transcriptomic, proteomic, and metabolomic levels, altering normal cellular metabolism and changing mitochondrial respiration (Wang et al. 2022).Disruption of normal mitochondrial function can result in an increase in reactive oxygen species (ROS) further exacerbating inflammation and increasing the likelihood of poor outcomes (Saleh et al. 2020).The "long COVID" phenomenon may be related to redox imbalance, which may be exacerbated by COVIDinduced changes in mitochondria (Paul et al. 2021;Singh et al. 2020) and, ultimately, fatigue related to metabolic impairment.Our results suggest that there is a delicate balance between metabolic demand changes and respiratory failure.One can easily speculate that reduction in available oxygen in concert with an increase in metabolic demand as the virus takes over cellular machinery to produce more viral particles can result in a point of critical failure.However, it is not intuitively obvious that greater O 2 carrying capacity results in metabolic collapse.Our model treats the relationship between oxygen carrying capacity and metabolic demand simplistically and we have not incorporated blood viscosity changes and their impact on cardiovascular function, particularly cardiac output as the key factor driving tissue perfusion.
Mitochondrial dysregulation due to COVID-19 results in pronounced impacts on blood coaguability (Hai-Han et al. 2020), handling of reactive oxygen species (ROS)increased by the cytokine storm associated with COVID-19 (Saleh et al. 2020), calcium homeostasis (Yang et al. 2021), iron homeostasis (Vlahakos et al. 2021;Abobaker 2020), as well as cellular metabolism (Henry et al. 2020;Booth et al. 2021).COVID-19 significantly impacts each of these aspects of mitochondrial function and this results in altered ability to respond to cytokine induced ROS changes, as well as reduced metabolic capacity.All of these alterations in mitochondria function contribute to changes in metabolic demand and it may be that, regardless of the O 2 concentration in the blood, the mitochondria are not able to utilize the available O 2 .In this study, we have shown that changes in metabolic reserve-particularly an impaired ability to meet metabolic demand in the context of respiratory functioncan result in collapse of respiration contributing to death.These changes are also most likely to be key factors in patients presenting with SH.
One way to test the impact of COVID-19 on mitochondrial function would be to assay mitochondria obtained from tissue biopsies of COVID-19 patients or through animal models.
Testing mitochondrial metabolism would be easier than using stress tests or cycle ergometry to determine metabolic load and ventilation/perfusion changes.Whole body tests would be problematic in COVID-19 patients and put them at greater risk for respiratory collapse.As long COVID has become better described, central nervous system (CNS) involvement and increased chronic inflammation are seen as sequelae that may continue to alter metabolism and mitochondrial function (Stefano et al. 2021).Further research is needed to determine if these effects are exacerbated by persistent metabolic impairment and whether symptoms like cognitive fog depend on CNS mitochondria and ROS handling problems.The relationship between changes in overall metabolic demand and cellular level metabolism have not yet been explored in COVID-19 patients.This is an important area for investigation because, while the metabolic demand required to pump more viscous blood (Choi et al. 2022) may selectively impact the cardiovascular system the most, metabolic demand may be increased systemically based on the diffuse organ involvement seen in these patients.
In addition to the model limitations we mentioned above, we realize that our model represents a very reduced number of the elements in the central pattern generator and pattern formation network for breathing control.The brainstem network includes hundreds of neurons that participate in each breath (Wang et al. 2014;Carroll and Ramirez 2013), and we have simplified this relatively complex circuit for the sake of rapid simulation time to test our hypotheses about SH.This heavily reductionist treatment of the brainstem network is an obvious limitation to simulation of the interacting populations of respiratory neurons and makes it difficult to interrogate the precise mechanisms by which respiratory collapse occurs in SH.Three examples include, (1) we have not included CO 2 sensing in our model due to the high diffusion rates of CO 2 when compared to O 2 in the lung (West 2008) and evidence showing that CO 2 is ≤ 35 mmHg in patients presenting with silent hypoxemia (SH) and minimal tachypnea (Chandra et al. 2020;Alamé et al. 2022); (2) we do not explicitly include rapidly adapting (RAR) or slowly adapting (SAR) lung mechanoreceptors in the model-lung volume is present in the model and reproduces inspiratory drive in much the same way that SARs do in vivo; (3) the lack of a specific mechanism for understanding the relationship between blood viscosity, number of red cells, oxygen carrying capacity, and changes in metabolic demand.
Previously, we have demonstrated that increasing extracellular [K + ] resulted in a progressive increase in respiratory rhythm that showed periodic, multi-periodic, quasi-periodic, and finally chaotic rhythmic patterns (Del Negro et al. 2002).As excitability increased, the disruption to eupneic breathing would result in impaired gas exchange in vivo.Thus, there is precedent for increasing excitability in the respiratory network resulting in a kind of "depolarization blockade" of normal breathing and a cessation of gas exchange that then results in a precipitous fall in P a O 2 .We described experiments related to this concept in our prior work (Diekman et al. 2017).Because we have previously shown these transitions are gradual and occur over a wide range of excitability changes, it makes sense to assume that there may be a more gradual progression of the "respiratory collapse," but we do not yet have clinical data showing how the collapse evolves to the point of need for ventilatory support.
Our model predicts changes in oxygen handling and metabolism in silent hypoxemia patients.As an obvious next step, we call for data to be collected on hematocrit in animal models of COVID-19 infection and COVID-19 patients, with the inclusion of metabolism and mitochondrial function tests.In addition, we believe the following measures may have untapped predictive value: minute ventilation, oxygen saturation, and breathing frequency.We speculate that some combination of these quantities, if measured on entry to the ER, could help predict the need for ventilator support in the subsequent 48 h.Finally, we emphasize that there is a need for incorporating oxygen handling dynamics into more sophisticated state-of-the-art respiratory control models, most of which currently focus on CO 2 and hypercapnea (Molkov et al. 2017).

A.1 Detailed model specification
Here, we provide the equations for the 7D-O2 model introduced in Diekman et al. (2017).Central pattern generator (CPG): A variety of models have been proposed for the central neural circuits generating breathing rhythms, ranging from group-pacemaker networks to individual pacemaker models, and beyond.Here, we adopt the original Butera-Rinzel-Smith (BRS) model (referred to as "model 1" in Butera et al. (1999a)) proposed as a mechanism for bursting pacemaker neurons in the preBötzinger complex.For simplicity, we represent the CPG with a single BRS unit.Thus, our CPG is described by a membrane potential V together with dynamical gating variables n (a delayed-rectifier potassium (I K ) channel activation) and h (persistent sodium (I NaP ) channel inactivation).We set two "instantaneous" gating variables p ∞ (I NaP activation) and m ∞ (fast sodium (I Na ) activation) to be equal to their voltagedependent asymptotic values.We set the I Na inactivation gate to be equal to (1 − n).The model also includes a leak current (I L ) and a tonic excitatory (I tonic ) current.In summary: it Motor pool activity The output of the CPG is the BRS cell's membrane potential (V ), which drives the respiratory muscles through synaptic activation of a motor unit (α): Here, r a = r d = 0.001 mM −1 ms −1 sets the rise and decay rate of the synaptic conductance.Also, [T ] represents the neurotransmitter concentration, with parameters T max = 1 mM, V T = 2 mV, and K p = 5 mV (Bard Ermentrout and Terman 2010).
Lung volume The output of the motor unit determines the rise and fall of lung volume (vol L ): Here vol 0 = 2 L is the volume of the unloaded lung, and parameters E 1 = 0.4 L and E 2 = 0.0025 ms −1 were chosen so that the lung expansion would remain in a physiologically reasonable range (West 2008).We note that while the low-frequency input of the envelope of CPG burst activity effectively drives changes in lung volume, highfrequency input (such as tonic spiking) does not drive the lung biomechanics effectively.This low-pass filter behavior of the respiratory musculature is analogous to tetanic muscle contraction that occurs in response to high-frequency stimulation of motor nerves (Kandel et al. 1991).
Lung oxygen At standard atmospheric pressure (760 mmHg), external air with 21% oxygen content will register a partial pressure of oxygen of P ext O 2 = 149.7 mmHg.As the lungs expand d dt [vol L ] > 0 , they draw in external air.Our model makes the simplifying assumption that this fresh air mixes instantaneously with the air already present in the lungs.
Therefore, the partial pressure of oxygen in the lung alveoli (P A O 2 ) increases at a rate given by the pressure difference between external and internal air, and by the lung volume.In contrast, during exhalation d dt [vol L ] ≤ 0 , no external air enters the lungs, so the mixing of air stops.During both contraction and expansion of the lung, oxygen moves between the lungs and the blood.The flux of oxygen from the lungs to the blood occurs at a rate determined by the time constant τ L B = 500 ms, and by the difference in partial pressure of O 2 between the lungs (P A O 2 ) and the arterial blood (P a O 2 ).The rate of change in P a O 2 is thus given by: Note the fluxes of oxygen from the blood to the tissues (J BT ) and from the lungs to the blood (J LB ) have units of moles of O 2 per millisecond.The term in the denominator converts from units of moles per millisecond (the rate of change of moles of O 2 in the blood) into units of millimeters of mercury per millisecond (rate of change of partial pressure of O 2 in the blood).To calculate the flux J LB , we use the ideal gas law PV = n RT , where n is the number of moles of O 2 , R = 62.364L mmHg K −1 mol −1 is the universal gas constant, and T = 310 K is temperature.The resulting flux depends on the difference in oxygen partial pressure between the lungs and the blood: Note that the term J BT accounts for both dissolved and hemoglobin-bound oxygen in the blood: The dependence of the hemoglobin saturation (SaO 2 ) on P a O 2 is given below, see Eqn. (24).Following Henry's law, we take the concentration of dissolved oxygen in the blood to be directly proportional to P a O 2 .The blood solubility coefficient, β O 2 = 0.03 ml O 2 × L blood −1 mmHg −1 for blood at 37 degrees C, is the constant of proportionality.The amount of dissolved O 2 at physiological partial pressures (P a O 2 ≈ 80 − 110 mmHg) is insufficient to satisfy the body's metabolic demand for oxygen.Therefore, most of the blood's stored oxygen is bound to hemoglobin (Hb).
Cooperative binding of oxygen to the four binding sites in each hemoglobin molecule leads to a sigmoidal hemoglobin saturation curve: Here, we take the phenomenological parameters K = 26 mmHg and c = 2.5 from Keener and Sneyd (2009).
Our model includes a parameter M in Eq. ( 23) to capture the rate of metabolic demand for oxygen from the tissues, in units of ms −1 .Equations ( 21) and ( 23 Chemosensation: Chemosensory feedback from peripheral chemoreceptors in the carotid bodies, carried to brainstem respiratory circuits via the carotid sinus nerve, close the control loop in our model.These receptors detect reductions in P a O 2 and drive the central rhythm generator, as described in more detail in Diekman et al. (2017).We model the nonlinear relationship between carotid chemosensory nerve fiber activity and P a O 2 as a sigmoid saturating function, with the firing rate low until P a O 2 is reduced below a threshold (normally about 100 mm Hg) and then steep firing rate increases as P a O 2 is reduced further (Hlastala and Berger 2001;West 2008).We capture this behavior in our model as a sigmoidal function connecting P a O 2 with the conductance representing external drive to the CPG (g tonic ): Here, φ = 0.3 nS, θ g = 85 mmHg, and σ g = 30 mmHg.This conductance closes the control loop in our respiratory control model, since I tonic = g tonic (V − E tonic ) is a term in the CPG voltage equation ( 7).We numerically integrated the preceding equations using a variable-order, variable-step stiff solver (ode15s in MATLAB).

A.2 Hemoglobin effects given the original model parameters
In this section, we illustrate the effects of changing the parameters controlling the hemoglobin binding curve and the total hemoglobin (hematocrit) in the original model, as opposed to the model adjusted by successive parameter shifts (see Sect. 3).Panel a shows the effect of increasing K , the hemoglobin binding affinity.Panel b shows the effect of varying [Hb], the hemoglobin concentration.The effects are qualitatively similar to the effects of analogous parameter changes in the silent hypoxemia model.

Fig. 2
Fig.2Dynamics of 7D-O2 model and hypothetical silent hypoxemia P a O 2 versus M curves.a Average blood oxygen P a O 2 as a function of metabolic demand (M) in the original 7D-O2 model.There is a plateau region at low M values for which normoxia (green shading) is maintained, and a collapse point at approximately M = 1.2 × 10 −5 ms −1 beyond which severe hypoxia occurs.In a model of silent hypoxemia, it seems clear that the plateau portion of the curve should shift lower (maroon arrow pointing down), but it is not as clear whether the collapse

Fig. 3 Fig. 5 aFig. 6
Fig. 3 Sensitivity of P a O 2 versus M curves to variation of model parameters.a Chemosensory sigmoid of g tonic as a function of P a O 2 with various parameter values for the maximum (φ), half-activation (θ g ), and slope (σ g ) of the sigmoid.Default settings from the original 7D-O2 model (φ = 0.3 nS, θ g = 85 mmHg, σ g = 30 mmHg) are shown in gray.See panel (b) for the definition of the color scale used for the other curves.b Average P a O 2 vs M curves for 25 different combinations of the chemosensory sigmoid parameters (φ = 0.24, 0.3, 0.36;

Fig. 2d .
Fig. 2d.Further increases in [Hb] yield collapse points with even lower M values, consistent with the hypothetical silent hypoxemia P a O 2 vs M curve shown in Fig. 2c.See Fig. 9b in Appendix A.2 for the effect of varying [Hb] with all other parameters set to their original 7D-O2 model values.Thus, we will consider the model with [Hb]=250 (the second curve from the left in Fig.5b) as a our working model for silent hypoxemia, and analyzed the model dynamics for simulations in the plateau region and in response to increases in metabolic demand.Figure6compares simulations of the silent hypoxemia model (blue traces) and the original 7D-O2 normoxia model (red traces) for different values of the metabolic demand parameter M, generated as follows.First, the corresponding model is simulated with M = 0.8 × 10 −5 ms −1 for two minutes of simulated time, to establish baseline initial conditions on the eupneic limit cycle.Then, the value of M is changed to the value shown on the horizontal axis (central column: panels b, e, h, k) and the simulation is run for another 10 min of simulated time.For the conditions shown in detail in the left column (M = 0.4 × 10 −5 ms −1 , panels a, d, g, j), this duration is sufficient to effectively remove tran- ) for the original 7D-O2 model and the SH model.Panels a (redrawn from Diekman et al. (2017)) and b show the phase line corresponding to the

Fig. 7
Fig. 7 Varying hemoglobin concentration shifts the collapse point in a reduced model obtained from fast-slow decomposition.a,b Rate of change of P a O 2 ( ḡ) of averaged slow subsystem (Eq.(3)) for M = .4× 10 −5 ms −1 (green), M = .8× 10 −5 ms −1 (blue) and M = 1.6 × 10 −5 ms −1 (magenta) for the models with normoxia (panel a) and silent hypoxemia (panel b) chemosensory parameters.Gray arrows in a indicate where P a O 2 lies in a range for which the fast subsystem sits in the quiescent state.c,d Heat map showing ḡ as a function of P a O 2 and M for the models with normoxia (panel c) and silent hypoxemia (panel d) chemosensory parameters.Black curves: fixed points ( ḡ = 0) of the averaged slow subsystem (3).e,f Fixed point ( ḡ = 0) curves for three levels of [Hb] (mmHg) for the normoxia (panel e) and silent hypoxemia (panel f) chemosensory parameters

Fig. 8
Fig. 8 Gray average P a O 2 versus M curves depict the mapping of the curves for 5 different values of hemoglobin concentration ([Hb] = 120, 180, 200, 250, and 300) to [Hb] = 150 by multiplying each curve by the ratio [Hb]/150 confirms this conjecture via simulations of the full model using several [Hb] levels.This figure replots the curves from Fig. 5b (colors from blue to red) together with rescaled versions of the figures (gray traces) obtained by multiplying each curve's abscissa by the ratio [Hb]/150, while keeping its original ordinate.As soon in the figure, the rescaled curves, plotted in gray, collapse onto approximately a common curve.
) include conversion factors ζ and η that depend on the concentration of hemoglobin, [Hb] = 150 gm L −1 , as well as the volume of blood, vol B = 5 L, respectively.The model assumes a molar oxygen volume of 22.4 L. We assume that each fully saturated hemoglobin molecule carries 1.36 ml of O 2 per gram: dt (P A O 2 ) = P ext O 2 − P A O 2 vol L where the notation [x] + indicates max(x, 0).Blood oxygen To represent the change in P a O 2 , we write