Theoretical open-loop model of respiratory mechanics in the extremely preterm infant

Non-invasive ventilation is increasingly used for respiratory support in preterm infants, and is associated with a lower risk of chronic lung disease. However, this mode is often not successful in the extremely preterm infant in part due to their markedly increased chest wall compliance that does not provide enough structure against which the forces of inhalation can generate sufficient pressure. To address the continued challenge of studying treatments in this fragile population, we developed a nonlinear lumped-parameter respiratory system mechanics model of the extremely preterm infant that incorporates nonlinear lung and chest wall compliances and lung volume parameters tuned to this population. In particular we developed a novel empirical representation of progressive volume loss based on compensatory alveolar pressure increase resulting from collapsed alveoli. The model demonstrates increased rate of volume loss related to high chest wall compliance, and simulates laryngeal braking for elevation of end-expiratory lung volume and constant positive airway pressure (CPAP). The model predicts that low chest wall compliance (chest stiffening) in addition to laryngeal braking and CPAP enhance breathing and delay lung volume loss. These results motivate future data collection strategies and investigation into treatments for chest wall stiffening.


Introduction
The extremely preterm infant, born at < 28 weeks gestation and often < 1000g, is at risk of developing chronic lung disease despite established treatments such as surfactant replacement therapy. Currently the survival rate of this group ranges from 94% at 27 weeks to as low as 33% at 23 weeks [1], with survivors living with varying degrees of morbidity. One risk factor for lung disease remains the trauma associated with traditional mechanical ventilation including endotracheal tube injury, high cyclic tidal volumes and pressures, and hyperoxia. Non-invasive methods of ventilation such as continuous positive airway pressure (CPAP) are being used with more frequency and have been successful with more mature infants but appear to fail in the extremely preterm infant [2][3][4]. One hypothesis for the failure of non-invasive ventilation a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 volumetric flow rate and rate of change represented as _ V ðtÞ [ml/s] and dV dt respectively. Air pressure P i within a specific volume i is defined as the difference between intra-airway pressure P int and pressure external to the body P ext , i.e. P i = P int,i − P ext,i . Since all pressures are relative to the same constant atmospheric pressure, all P ext = 0 and all intra-airway pressures P i = P int,i . The pressure P ij = P i − P j refers to the transmural pressure across a compliant boundary separating volumes i and j.

State equations
Each non-rigid compartment has an associated compliance C i [ml/cm H 2 O], describing the change in compartmental volume V i given a change in transmural pressure P ij across its boundary with compartment j: The nature of C i does not change explicitly with time but instead is implicitly determined by  the relationship between volume and pressure. This can be reformulated in terms of dynamic changes of state: Bidirectional airflow through the trachea, bronchi, bronchioles, and to and from the lungs results from contraction and relaxation of the diaphragm generating a pressure difference. Airflow is opposed by the resistance of the airways as functions of their radaii or tissue properties. This relationship is described by the flow-pressure analog of Ohm's law [19], where R i [cm H 2 OÁs/ml] is the resistance to airflow prior to compartment i. If a compartment includes inertial effects, the pressure gradient is also a function of the acceleration of flow, where I is the inertance. Inertial effects are considered for the newborn upper rigid airway because of its smaller radius, but neglected for the rest of the model tissues [17]. The pressures P ij across each compliant compartment include transmural pressure between the compliant airways and the pleural space P tm = P c − P pl , lung elastic recoil P el = P A − P T , lung viscoelastic component P ve = P T − P pl , and chest wall elastic recoil P cw = P pl − P mus . Summing pressures over each of three loops according to Kirchhoff's mesh rule gives a system of time-varying algebraic equations: Using the additional relationships obtained from applying Eqs (1) and (2), the system of loop equations can be rewritten as Rearranging Eq (3) and using Kirchhoff's current law along with Eqs (1) and (2) produces the consolidated set of model differential equations: Conservation laws also maintain that V = V cw = V A + V c , in other words the total system volume equals the chest wall volume, which is the sum of the alveolar and compressible airway volumes. Pressure-volume relationships and compliances C i will be further described below.

Nonlinear resistance constitutive relations
The airways begin with an upper rigid segment characterized by an inertance I u and a nonlinear Rohrer resistance R u [17,20] that increases with airflow: The constants R u,m and K u represent laminar and turbulent flow components. A middle collapsible portion is modeled as a cylinder with constant length having nonlinear resistance R c that depends inversely on the 4th power of the radius according to Poiseuille's law. Therefore R c is formulated as [12,21]: where R c equals its minimum value K c when V c = V c,max , an estimate of dead space. An inverse relationship between resistance in the smaller peripheral airways R s and lung volume V A reflects high resistance at low or near-zero volumes [21]. To avoid R s ! 1 as V A ! 0 [22] from a strict exponential decay model, we adopt the formulation used by both Liu et al and Athanasiades et al [12,14], a decaying exponential function of relative lung volume with finite R s at V A = 0: where K s < 0. This parameterization gives that R s % R s,m when V A = TLC, and R s = R s,d + R s,m when V A = RV (residual volume).

Nonlinear compliance constitutive relations
The volumes V c , V A , and V cw representing physiological compartments are assumed to have nonlinear compliance which are modeled implicitly with a pressure-volume curve or explicitly by C i ¼ dV i dP i . The compliance curve for the collapsible airway volume V c as a function of P tm represents data depicted in [21] following a sigmoidal function [17,22] Maximal compliance occurs at the middle of the sigmoid c c , with d c characterizing the slope of the sigmoid. In newborns and infants an exponential-like chest-wall compliance curve is observed [23,24] but with compliance being near infinite for P cw > 0. We chose to model the static compliance of the chest wall as a "softplus" function of the form f(x) = ln(1 + e x ), the smooth approximation of the rectifier activation function f(x) = max(0, x). Accounting for translations and scaling, this is represented by The asymptotic volume at large negative pressure is thus assumed to equal RV. The "transition point" where the softplus function slope has the greatest rate of change from horizontal to affine occurs at P cw = . The chest wall relaxation volume V 0 = V cw | P cw = 0 is set using an estimate from literature at 25% of VC (vital capacity) [23,24]. From this parameterization, b w = (V 0 − RV)/(ln 2). The single degree of freedom d w then characterizes the slope of the chest wall compliance curve and is adjusted to produce a range of dynamic compliance values. The volume of the lung compartment V A is modeled as the product of distention of lung units V el (P el ) and fraction of recruited alveoli F rec (P el ) [25,26]. To obtain V A % RV near P el = 0, lung volume is given as Alveolar compliance C A as used in the system of differential Eq (4) is found with symbolic computation as dV A dP el . The first term V el represents the volume due to aggregate elasticity of the lung unit structure, which is modeled here as a saturated exponential [26][27][28] where k characterizes the lung stiffness. This representation has been found to suffice in cases of a healthy or surfactant-treated lung. The second term of the lung compliance F rec represents the contribution of recruitment and derecruitment of alveoli to compliance, which has been modeled previously as dependent on both time and pressure [25,26,29]. It can be represented by a sigmoid which resembles the probability density function of a Gaussian distribution describing aggregate opening or closing pressures of individual alveolar sacs or ducts [30]. We adopt the formulation of Hamlington et al [26]: It follows that β is the baseline fraction of lung recruited at P el = 0, γ represents the maximum recruitable fraction of lung, c F is mean opening pressure at which recruitment is maximum, and d F describes the transition to full recruitment capturing the heterogeneity of the lung. Parameterization of F rec is based on the state of health being modeled and can change breath-to-breath depending on conditions. For example, an increase in stiffness resulting from derecruitment may manifest as higher mean opening pressure c F and move the V A curve to the right. Likewise a lower maximum recruitable fraction γ would flatten the V A curve. Both scenarios indicate a lower compliance and greater pressure required to increase the lung volume in the region of operating pressure. In certain pathological situations such as ARDS, a sigmoidal representation of V A (P el ) with a low compliance region at low P el [28,[31][32][33][34] could be captured in the parameterization of F rec .
The viscoelastic properties of pulmonary tissue are represented with a linear Kelvin-Voigt model consisting of scalar compliance C ve and resistance R ve , which contributes a viscoelastic pressure component P ve in series with lung elastic recoil P el , see Fig 1. The sum of these two pressures is dynamic pressure P l,dyn which also equals P pl − P A .

Respiratory muscle driving pressure
The pressure P mus describes the effective action of the respiratory muscles driving the model dynamics with P mus negative in the outward direction. We used a sinusoidal function to describe tidal breathing, with maximum equaling zero at end-expiration: where A mus is the amplitude of the cosine wave and f = RR/60 is the frequency. The wave generates a negative pressure with total magnitude 2A mus outward from the body. Though simple, the sinusoidal function can admit time-varying frequency, show dynamics over multiple breaths, is used in artificial ventilation, has compact support on the closed interval [0, T], and has been used in previous modeling studies (see eg. [14]). More sophisticated functions [17,35,36] can model inhalation and exhalation with different durations or qualitative forms, however the breath-to-breath dynamics displayed in this study can be captured sufficiently with the sinusoidal function.

Progressive volume loss
The complete mechanism of interaction between inefficient inhalation resulting from high chest wall compliance and the progressive nature of lung volume loss and respiratory distress is not fully understood. Clinical X-ray evidence of delayed atelectasis and subsequent acute respiratory distress in otherwise healthy lungs may suggest a process by which a lack of full recruitment during a given breath lowers lung capacity and compliance for the following breath, and continues to an unrecoverable level in the absence of neural modulation or compensatory mechanisms such as sighing. As a first attempt at modeling progressive volume loss, we empirically describe the breath-to-breath evolution of F rec (Eq (12)) as lung recruitment pressure parameters c F and d F increase with PIP and maximum recruitable fraction γ decreases with EILV. The lung compliance curve shifts slightly with each breath via changes in mean threshold opening pressure based on number of collapsed alveoli. A volume loss associated with derecruited alveoli necessitates an increase in expanded volume of recruited alveoli relative to the radius cubed, with an increased distending pressure proportional to the change in radius. This is illustrated in [37] using a simple example of expansion of 3 alveoli that double in volume with a 25% increase in radius; if 1 alveolus closes, the other two radaii must now increase by 35% to achieve the same overall volume change and the required distending pressure increases proportionally. This proportion applied to c F and d F shifts the compliance curve to the right. In this way the compliance decreases approximately proportional to the amount of derecruitment [38]. If tidal breathing begins on the steepest part of the lung compliance curve, compliance decreases monotonically until eventually tidal breathing occurs on the low compliance tail on the left part of the curve and V T % 0. Tidal breathing may begin at a higher position towards the flatter upper part of the curve, in which case compliance will increase slightly with this modification but will again eventually decrease in the manner described above.
Assuming constant amplitude of the sinusoidal muscle pressure pressure function and no stochasticity, the maximum recruitable fraction of alveoli is achieved at end-inspiration (EI) during steady-state oscillatory breathing and additional fraction will not be recruited under a pressure of this same amplitude in subsequent breaths. The value for γ for subsequent breaths is then dependent on F rec | EI and the percentage of alveoli assumed to be permanently collapsed / no longer recruitable, represented by the calculation g next ¼ g current _ ð1 À ð% permanent Þ Á ðF closed ÞÞ where F closed is fraction closed. If all unrecruited alveoli remain as such, then F rec | EI becomes the new γ for the next breath; likewise, if all alveoli remain recruitable, γ = 1 for the duration of the simulation. Note that even for γ = 1, F rec < 1 for all P el thus causing small changes in c F and d F and shifts in the F rec (P el ) curve regardless of the % of alveoli permanently closed that still lead to progressive volume loss. The rate at which volume loss progresses depends on where on the F rec curve tidal breathing occurs, and thus both the curve's intrinsic characterizing parameters and extrinsic system variables.

Parameterization
The lung curve was parameterized to obtain an approximate dynamic lung compliance C A of 2.3 ml/cm H 2 O [8,39,40] calculated as the slope (V A | EI − V A | EE )/(P el | EI − P el | EE ) during normal breathing with no interventions. In particular, k was tuned to produce a curve V el between RV and TLC with the calculated slope, and the parameters of F rec produced a curve that is % 1 for the whole range of normal breathing to represent a nearly fully recruited lung. High C w for a typical preterm infant was targeted at 8.5 ml/cm H 2 O [8,39] and low C w about equal to lung compliance. The parameter d w characterized the approximate dynamic chest wall compliance, which was calculated as the slope (V cw | EI − V cw | EE )/(P cw | EI − P cw | EE ). Parameter values for R u , R c , R s , and V c were estimated from previously published studies [12,14,41]. The viscoelastic parameters C ve and R ve were manually tuned to obtain idealized tidal volume and end-expiratory lung volume rather than the magnitude of the hysteresis.
FRC is the volume at the resting position of the respiratory system i.e. where P resp = P el + P cw = 0. The naturally high compliance of the healthy full-term and especially preterm infant (with even steeper V cw (P cw )) lowers P resp and decreases FRC to about 20% of vital capacity (VC), compared to at about 35-40% of VC in the adult [42]. A nominal value for FRC for a given set of static compliance curves is obtained by first computing volumes using a vector of physiological pressures [-20. . .40] cm H 2 O. Lung and chest recoil pressure vectors are then added in the P direction to obtain P resp , and the index where P resp = 0 is used to determine FRC using either V cw or V A . The lung, chest wall, and respiratory PV compliance curves for both high and low C w created from Eqs (9) and (10) are given in Fig 2. The value for P el | FRC is then Fig 2. Lung, chest wall, and total respiratory system compliance curves for high C w (left) and low C w (right). Curves are described by Eqs (9) and (10) and parameterized using the procedures described in Parameterization. Tidal breathing loops with normal R u (grey) and increased R u (black) are superimposed for each condition over the lung compliance curve and larger in each inset to display hysteresis. set as the initial condition for solving dP el /dt. Note that the lung curve is identical between scenarios so the decreased slope in the low C w scenario with the same ν and V 0 raises FRC and thus EELV. Decreased lung compliance (flatter V A (P el )) resulting from injury, disease, or progressive volume loss further reduces FRC and EELV. In this model we consider chest wall compliance to be either high or low and unchanging for the duration of a simulation, but lung compliance changes depending on breathing conditions. Table 2 gives values and formulas / sources for parameters that remain unchanged between simulations. These values as well as the FRC, respiratory pressure amplitude, chest wall compliance, and upper airway resistance parameters in Table 3 that vary between simulation conditions were manually tuned to best obtain the reported aggregate parameters and state outputs as shown in Table 4. As an example, dynamic lung compliance C A is not an explicit input into the model, but was determined as described above. For ease of computation and to match the target demographic, we assumed the simulated subject weighed 1 kg.

Computational procedures
All simulations proceeded with an initial respiration rate of 60 breaths/min (f = 1), initial minute ventilation _ V E ¼ 360, and initial tidal volume V T = 6 ml, with the expectation that tidal volume changes with changes in dynamic lung compliance. The motivation for ths choice was twofold: One, this is consistent with a physiological requirement of constant _ V E regardless of chest or lung compliance; and two, this allowed for comparison of simulation results originating from similar starting points. Distinct values for A mus were prescribed for each simulation to achieve the initial _ V E ¼ 360, see Table 3. Simulated conditions were chosen to demonstrate the model dynamics with high and low chest wall compliance, under two interventions and two states of permanent alveolar closure. An infant often exhibits compensatory mechanisms such as laryngeal braking (grunting) and increased activity of diaphragm and intercostal muscles [18] to increase end-expiratory pressure in order to keep EELV above the volume at which alveolar units start to collapse during expiration. Laryngeal braking is simulated with a 10-fold increase in expiratory upper airway resistance R u . CPAP is simulated with an increase of P ao from 0 to 5 triggered at F rec,EI = 0.9, 0.95, 0.97 characterizing volume losses of 10%, 5%, and 3%. Simulations also include assumptions of either no permanently closed alveoli, such that γ = 1 for all time, or10% permanently closed alveoli per breath, such that g next ¼ g current _ ð1 À 0:1ðF closed ÞÞ. Each simulation was ave is a moving average of the previous 60 tidal volumes (% 1 minute of breathing); variable f including a single 20 second apneic event.
The system of differential Eq (4), together with the constitutive relations (5-13), were solved using MATLAB R2016b (MathWorks, Natick, MA) with the differential equations solver ode15s (see S1 Code). Initial conditions were set at physiological values as given in Table 5. The equations were solved for each new breath using the end conditions from the previous breath as initial conditions. Parameter values as discussed earlier are given in Tables 2  and 3. The steady-state stability of the model was analyzed under constant non-oscillatory muscle pressure by examining the eigenvalues of the Jacobian at the nominal parameter set and varying parameters by multiples of 2 and 10. Results of this analysis are found in S1 Appendix.

Results
Parameterized static compliance curves for V cw (P cw ) and V A (P el ) are shown in Fig 2 for high C w (left) and low C w (right). The hysteretic tidal breathing loops are superimposed on the curve V A (P el ) for normal R u in black and increased R u in grey. Hysteresis is caused in the model by the viscoelastic parameters C ve and R ve , which were tuned to maintain appropriately valued lung volume outputs. Fig 3 shows the impact of high vs low C w and normal vs. high R u on the five states P A , P l,dyn , P pl , V A , and _ V . Increased R u increases P A almost threefold, but C w has very little impact. However, decreased C w increases P l,dyn significantly, effectively raising it higher on the lung PV curve. Increasing R u even higher increases P l,dyn but there is no difference with respect to C w . The opposite appears to occur with P pl dynamics, in that decreasing C w makes P pl more negative ("increasing" the magnitude of the pressure) and increasing R u strengthens that effect. Low C w and subsequently high R u increase V A , mimicking the effect for P l,dyn . High R u shifts _ V by 5 ml/s, with airflow more restricted during expiration. Tabulated magnitudes of the steady states are gives in Table 4. These results compare favorably to the record reported in Abbasi et al [50] in which esophageal (pleural) pressure, airflow, and tidal volume were approximately -2 to -6 cm H 2 O, -30 to 30 ml/s, and 8 ml, respectively. Table 6 presents the 14 simulations and their time to failure, defined for this study as 90% volume loss. Dynamics were comparable between simulations with the major difference being the timing, therefore only representative or significant results are presented in figures. Our model consistently indicates a faster loss of end-expiratory lung volume in all simulations with high C w compared to the same with low C w . Variable f did not significantly change TTF except in the case of CPAP administered at 3% loss (S14), with TTF shortened by almost 2 hours. Adding a single 20 second apnea shortened the TTF by 1-4 minutes in the shortest simulations but by over an hour under high C w and increased R u (S8). The breath-to-breath change in EELV and V T under high and low C w conditions with no interventions are given in Fig 4 (Simulations 1 and 3). The high C w simulation reaches accelerated loss of volume and eventually failure at 0.3 hours, much more quickly than the low C w at 2.5 hours. This depicts a possible scenario in which lung volume loss and failure may appear to onset suddenly after a long period of apparent steady conditions. Fig 5 shows changes in dynamic lung compliance and tidal volume with high and low C w without changes in γ, then adding CPAP to the high C w condition at three different levels (c.f. Table 6, simulations 1, 3,11,[13][14]. CPAP was simulated by an increase in mouth pressure P ao to 5 cm H 2 O when F rec,max < 0.9, which happened when the lung volumes were already decreasing quickly towards failure. However, CPAP triggering at F rec,max < 0.95 and 0.97 gained 3 and 9 hours of time, respectively. Note that regardless of timing, the administration of CPAP is correlated with reduced tidal volume (see also [51]). Increasing P ao moves the resulting PV loop higher up on the lung compliance curve but does not change the nature of the curve, thus eventually the influence of high C w on dynamics induces the same lung volume loss without other mitigating actions.

Discussion
In summary, we have developed a lumped-parameter respiratory mechanics model tuned with parameters specific to the extremely preterm infant weighing 1 kg. The model includes a novel representation of derecruitment based on alveolar pressure and volume expansion compensating for collapsed alveoli. Model simulations suggest conditions under which volume loss may result more quickly from higher vs lower chest wall compliance in the preterm infant, indicating the plausibility of dynamics underlying the symptoms observed clinically. Given the fragile nature of this population, it is extremely difficult to obtain non-pathological parameter or state output values for a healthy or surfactant-treated infant during spontaneous breathing, and even more so to obtain time series for model validation and eventual parameter estimation. The much earlier study by Abbasi and Bhutani [50] and a later one by Pandit et al [40] gave the best insight into the respiratory dynamics of an extremely preterm infant, making these the standard against which our results were qualitatively validated. We therefore claim that this effort is a "proof of concept" that will be further explored in future investigations using pressure and airflow time series data in a parameter estimation / optimization procedure to characterize parameter values specific to a particular patient dataset. Additional model modifications will allow for hypothesis generation for future data testing and data collection. As mentioned briefly in Nonlinear compliance constitutive relations, recruitment/decrecruitment may have a time component [29,52], in that the time it takes for an airway or alveolus to open may be a function of how far away its pressure is from its critical opening pressure. Earlier studies have developed models that incorporate opening and closing pressures for individual alveoli, contributing to the aggregate difference in inflation and deflation limbs of the hysteretic PV curve [53,54]. These previous studies considered recruitment resulting from one or two hyperinflations but not long-term derecruitment. In our model breath-to-breath derecruitment is manifested as the change of the lung compliance curve during normal spontaneous breathing as described in Progressive volume loss, and the hysteresis found in the tidal breathing loop is accounted for by the viscoelastic component of the system of differential equations. It is clear from Table 6 that time to failure shortens if an assumption is made about a non-zero percentage of alveoli permanently closing and being unavailable for recruitment. As a topic for further study, the pulmonary tissue may be modeled by more complex Voigt-Maxwell models within the "electrical analog" model or other non-electrical analog representations (e.g. those described in [29]). While such a modification may affect the overall trends in observed states such as EELV, the differential impact between high and low C w would be expected to remain.
The noninvasive ventilatory intervention CPAP shifts the tidal volume loop to a higher position on the lung compliance curve, operating with a higher EELV and end-expiratory lung elastic recoil. Our model suggests that the timing of administration of CPAP and the permanent closure or injury state of alveoli may impact its effectiveness. In our first simulation with simulated CPAP triggered at 10% volume loss, the recruited volume fraction does not recover fully to 1 and the use of CPAP only gains about a half hour of breathing before failure. However, CPAP starting at 5% and 3% loss gained 3 and 9 hours of time, respectively. This magnitude of loss may not be symptomatic at this point but would benefit from pressure support to avoid the quick descent to failure. These results are reported for the case with all fully recruitable alveoli. In the case of 10% permanent collapse of closed alveoli at each breath and subsequent breath-by-breath decrease in γ, the function F rec can never reach 1 (full recruitment) for the duration of breathing, tidal breathing occurs on a lower lung compliance curve, and CPAP cannot recover the full volume loss in subsequent breaths. Results in Table 6 indicate that time to failure is 10% faster with the permanent collapse. These simulated loss and collapse percentages were arbitrarily chosen to demonstrate the capabilities of the model and possible influences on breathing dynamics, but more investigation into actual loss values would add to the model's usefulness. Starting from a lower C w appears to be the optimal condition presented here as hypothesized.
Prolonged shallow breathing has been associated with increased surface tension and decreased surface area that further hinders breathing [55]. We safely assume in our model that derecruitment is a continuous process that will eventually induce loss of lung volume if left uncompensated [56,57]. In a healthy lung in the absence of fatigue, permanent alveolar collapse (due to injury or disease), and/or high chest wall compliance, this process is on a much longer time scale than the natural compensation mechanisms that compensate for and recoup volume loss (such as grunting in the infant). One such mechanism is spontaneous deep breathing, or "sighing", which may help prevent atelectasis [58][59][60] by re-opening air spaces that collapse naturally under tidal breathing [61] via increased pressure and surfactant activation and possibly affect neurorespiratory control. Sighing occurs more frequently and at relatively larger magnitude in the infant vs adults [62]. A natural extension of our model would be incorporating the restorative actions of sighing and testing the hypothesis that spontaneous deep breaths mitigate or reverse volume loss.
Several features of the physiology of preterm infants are not currently addressed in this model but should be considered in future model enhancements for further investigations. Preterm infants commonly exhibit diaphragm weakness and dysfunction and paradoxical breathing. While a sinusoidal waveform is used in the clinic under some mechanical ventilation protocols, the sinusoidal pressure function used here is an elementary representation for spontaneous breathing and does not capture dynamics related to diaphragm dysfunction or possible expiratory flow limitation. Modifications reflecting such dynamics may include adjustments to the pressure amplitude, varying fractions of time spent in inspiration vs expiration, and the use of a model that combines functional forms such as polynomials or exponentials (see e.g. [35]). Components that differentiate between abdominal and rib cage movements (see e.g. [63]) may model the paradoxical chest movement.
Another limitation of this model is the absence of any feedback mechanisms compensating for loss of volume. More sophisticated models of central pattern generators have been developed in conjunction with simple lung mechanics [64,65] that could potentially be incorporated with ours. A chemoreflex model, see for example [65,66], may also augment our model. Despite these limitation, we expect that the timing of dynamics of individual simulations may change with model enhancements but that time to failure would still be extended under low chest wall compliance conditions as observed in this study.

Conclusion
Respiratory mechanics models have been investigated for several years and many formulations exist; the challenge to be appreciated is the customization to the preterm infant with significantly different physiological features than adults and even term infants. Hence future model modifications must always keep this at the forefront of any investigation. The lumped-parameter respiratory mechanics model developed in this study will be used in future studies with data currently being collected in the NICU to estimate patient-specific parameters, which may shed light on factors influencing volume loss dynamics. This process may help generate hypotheses about predicting volume loss and recovery to motivate future data collection strategies. Our hope is that these investigations lead to a chest-stiffening treatment that can target an infant's specific physiological characteristics and prevent volume loss in this vulnerable population.
Supporting information S1 Appendix. Model stability analysis. The inherent stability of the model was analyzed under constant non-oscillatory muscle pressure by examining the eigenvalues of the Jacobian at the nominal parameter set and varying parameters by integer multiples. (PDF) S1 Code. Minimal representative code. Attached MATLAB code runs Simulation 3 (high C w , normal R u ) with constant frequency for 6 periods. (PDF)