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-4]. One hypothesis for the failure of non-invasive ventilation and the need for increasing invasive respiratory support is the markedly increased compliance (floppiness) of the chest wall in the extremely preterm infant resulting from ribcage undermineralization common at the start of the third trimester [5-7]. In the preterm infant, chest wall compliance can be up to five times lung tissue compliance [8].

When the chest wall is not sufficiently rigid, the negative pressure within the pleural space between the lung and chest wall generated from diaphragm contraction is diminished [9]. In many cases this leads to progressive lung collapse (atelectasis) with each breath as the forces needed to open airspaces after each exhalation become insurmountable [10], leading to decreasing lung compliance and functional residual capacity (FRC) [11]. This progression of events is observed clinically in X-rays and by symptoms of respiratory distress such as chest retractions and rapid breathing. The clinical result is progressively reduced tidal volumes and end-expiratory lung volume (EELV) as the forces needed to open airspaces after exhalation are insufficient. Non-invasive ventilation has been observed to be become ineffective under these conditions, necessitating placement of an endotracheal tube and positive pressure mechanical ventilation and markedly increasing the risk of lung damage.

Despite this being repeatedly observed clinically, there remains little quantification of the impact of variable nonlinear chest wall compliance on tidal breathing dynamics, and even fewer computational modeling efforts supporting these observations of progressive volume loss. Most computational models of breathing address the extremes of lung capacity such as a forced vital capacity maneuver, study a static, excised, or injured lung, or use an animal model [12-16]. Existing computer models of tidal breathing have not fully accounted for the physiology particular to premature infants and thus have limited apan class="Gene">pplicability. Often, methods of providing ventilator support have been developed in adults and children, then refined and scaled for newborns and premature infants, limiting innovation aimed specifically at this vulnerable population.

In this work, we have developed a nonlinear computational model of respiratory mechanics parameterized for the extremely preterm infant that demonstrates differential volume loss under high vs low chest wall compliance conditions. We adapt a model first presented by Athanasiades et al [14] and modified for newborn lambs by LeRolle et al [17]. In the latter, differences such as smaller diameter airways, higher respiratory rates, higher lung resistance, and higher chest wall compliances were considered, however many of the critical physiological nonlinearities contributing to long-term dynamics were not included. The present model is built upon the nonlinear compliance curves describing pressure-volume relationships specific to preterm infants [18]. Dynamic alterations of compliance curves based on breath-to-breath end-inspiratory lung volume (EILV) and peak inspiratory pressure (PIP) are shown to influence tidal volume and EELV, thus simulating progressive lung volume loss. We also demonstrate the effect of two simulated interventions that raise alveolar pressure and lung elastic recoil: CPAP, which raises the pressure at the mouth; and laryngeal braking (grunting), which increases upper airway resistance during expiration.

Mathematical model

The lumped-parameter respiratory mechanics model describes dynamic volumes and pressures in the airways, lungs, chest wall, and intrapleural space between lungs and chest. A signal that represents diaphragm pressure generated during spontaneous breathing drives the model. A compartment is assumed to display aggregate behavior, e.g. the alveolar compartment represents the collective dynamics of the alveoli as a whole. The model is designed using the volume-pressure analog of an electrical circuit, see Fig 1. As such, relevant states are in terms of pressure P(t) [cm H2O] and volume V(t) [ml] in and between air compartments, with volumetric flow rate and rate of change represented as [ml/s] and respectively. Air pressure P within a specific volume i is defined as the difference between intra-airway pressure P and pressure external to the body P, i.e. P = P − P. Since all pressures are relative to the same constant atmospheric pressure, all P = 0 and all intra-airway pressures P = P. The pressure P = P − P refers to the transmural pressure across a compliant boundary separating volumes i and j.

Fig 1

Lumped-parameter respiratory mechanics model, in both volume-pressure (panel A) and electrical (panel B) system analogs.

Each non-rigid compartment has a volume V (black), pressure P, (black) and associated compliance C (green, for emphasis) that is a function of the transmural pressures (purple) across the compartment boundaries. Air flows (red) across resistances R and inertance I (blue) are positive in the direction of the arrows. Circular yellow arrows indication direction of loop summations in Eq (3). Subscripts: airway opening ao, upper u, collapsible c, small peripheral s, alveolar A, viscoelastic ve, lung elastic el, transmural tm, pleural pl, chest wall cw, muscle mus.

State equations

Each non-rigid compartment has an associated compliance C [ml/cm H2O], describing the change in compartmental volume V given a change in transmural pressure P across its boundary with compartment j: The nature of C 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 [cm H2O⋅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 across each compliant compartment include transmural pressure between the compliant airways and the pleural space P = P − P, lung elastic recoil P = P − P, lung viscoelastic component P = P − P, and chest wall elastic recoil P = P − P. 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 = V + V, 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 will be further described below.

Nonlinear resistance constitutive relations

The airways begin with an upper rigid segment characterized by an inertance I and a nonlinear Rohrer resistance R [17, 20] that increases with airflow: The constants R and K represent laminar and turbulent flow components.

A middle collapsible portion is modeled as a cylinder with constant length having nonlinear resistance R that depends inversely on the 4th power of the radius according to Poiseuille’s law. Therefore R is formulated as [12, 21]: where R equals its minimum value K when V = V, an estimate of dead space.

An inverse relationship between resistance in the smaller peripheral airways R and lung volume V reflects high resistance at low or near-zero volumes [21]. To avoid R → ∞ as V → 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 at V = 0: where K < 0. This parameterization gives that R ≈ R when V = TLC, and R = R + R when V = RV (residual volume).

Nonlinear compliance constitutive relations

The volumes V, V, and V representing physiological compartments are assumed to have nonlinear compliance which are modeled implicitly with a pressure-volume curve or explicitly by .

The compliance curve for the collapsible airway volume V as a function of P represents data depicted in [21] following a sigmoidal function [17, 22] Maximal compliance occurs at the middle of the sigmoid c, with d 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 > 0. We chose to model the static compliance of the chest wall as a “softplus” function of the form f(x) = ln(1 + e), 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 = . The chest wall relaxation volume V0 = V| is set using an estimate from literature at 25% of VC (vital capacity) [23, 24]. From this parameterization, b = (V0 − RV)/(ln 2). The single degree of freedom d 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 is modeled as the product of distention of lung units V(P) and fraction of recruited alveoli F(P) [25, 26]. To obtain V ≈ RV near P = 0, lung volume is given as Alveolar compliance C as used in the system of differential Eq (4) is found with symbolic computation as .

The first term V represents the volume due to aggregate elasticity of the lung unit structure, which is modeled here as a saturated exponential [26-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 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 = 0, γ represents the maximum recruitable fraction of lung, c is mean opening pressure at which recruitment is maximum, and d describes the transition to full recruitment capturing the heterogeneity of the lung. Parameterization of F 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 and move the V curve to the right. Likewise a lower maximum recruitable fraction γ would flatten the V 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(P) with a low compliance region at low P [28, 31–34] could be captured in the parameterization of F.

The viscoelastic properties of pulmonary tissue are represented with a linear Kelvin-Voigt model consisting of scalar compliance C and resistance R, which contributes a viscoelastic pressure component P in series with lung elastic recoil P, see Fig 1. The sum of these two pressures is dynamic pressure P which also equals P − P.

Respiratory muscle driving pressure

The pressure P describes the effective action of the respiratory muscles driving the model dynamics with P negative in the outward direction. We used a sinusoidal function to describe tidal breathing, with maximum equaling zero at end-expiration: where A is the amplitude of the cosine wave and f = RR/60 is the frequency. The wave generates a negative pressure with total magnitude 2A 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 (Eq (12)) as lung recruitment pressure parameters c and d 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 and d 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 ≈ 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| and the percentage of alveoli assumed to be permanently collapsed / no longer recruitable, represented by the calculation where F is fraction closed. If all unrecruited alveoli remain as such, then F| 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 < 1 for all P thus causing small changes in c and d and shifts in the F(P) 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 curve tidal breathing occurs, and thus both the curve’s intrinsic characterizing parameters and extrinsic system variables.

Simulation conditions

Parameterization

The lung curve was parameterized to obtain an approximate dynamic lung compliance C of 2.3 ml/cm H2O [8, 39, 40] calculated as the slope (V| − V|)/(P| − P|) during normal breathing with no interventions. In particular, k was tuned to produce a curve V between RV and TLC with the calculated slope, and the parameters of F produced a curve that is ≈ 1 for the whole range of normal breathing to represent a nearly fully recruited lung. High C for a typical preterm infant was targeted at 8.5 ml/cm H2O [8, 39] and low C about equal to lung compliance. The parameter d characterized the approximate dynamic chest wall compliance, which was calculated as the slope (V| − V|)/(P| − P|). Parameter values for R, R, R, and V were estimated from previously published studies [12, 14, 41]. The viscoelastic parameters C and R 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 = P + P = 0. The naturally high compliance of the healthy full-term and especially preterm infant (with even steeper V(P)) lowers P 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 H2O. Lung and chest recoil pressure vectors are then added in the P direction to obtain P, and the index where P = 0 is used to determine FRC using either V or V. The lung, chest wall, and respiratory PV compliance curves for both high and low C created from Eqs (9) and (10) are given in Fig 2. The value for P| is then set as the initial condition for solving dP/dt. Note that the lung curve is identical between scenarios so the decreased slope in the low C scenario with the same ν and V0 raises FRC and thus EELV. Decreased lung compliance (flatter V(P)) 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.

Fig 2

Lung, chest wall, and total respiratory system compliance curves for high C (left) and low C (right).

Curves are described by Eqs (9) and (10) and parameterized using the procedures described in Parameterization. Tidal breathing loops with normal R (grey) and increased R (black) are superimposed for each condition over the lung compliance curve and larger in each inset to display hysteresis.

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 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.

Table 2

Tuned steady-state and dynamic simulation parameters that remained unchanged during simulations.

Parameter	Value	Formula	References
TLC [ml]	63	—	[18, 23]
RV [ml]	23	—	[18]
VC[ml]	40	TLC-RV	[18, 23]
RR [br/min]	60	—	[23]
f [br/s]	1	RR/60	—
T [s]	1	1/f	—
ν	0.25	—	[23, 24]
V0 [ml]	35	ν⋅VC+RV	—
β	0.01	estimated	[26]
γ	1	estimated	[26]
α	-0.76	(1+e^c_F/d_F)β-γe^c_F/d_F	[26]
k [1/cm H2O]	0.07	estimated	[26, 28]
cF [cm H2O]	0.1	estimated	[26]
dF [cm H2O]	0.4	estimated	[26]
aw [ml]	23	RV	[23, 24]
bw [ml]	17.3	(V0 − RV)/ln2	—
cw [cm H2O]	0	estimated	—
ac [ml]	0	—	[12]
bc [ml]	2.5	Vc,max	[43]
cc [cm H2O]	4.4	estimated from adult	[12]
dc [cm H2O]	4.4	estimated from adult	[12]
Kc	0.1	estimated from adult	[14]
Vc,max [ml]	2.5	estimated as dead space	[23, 43]
Rs,m [cm H2O s/L]	12	—	[41, 44]
Rs,d [cm H2O s/L]	20	estimated from adult	[14]
Ks	-15	estimated from adult	[14]
Iu [cm H2O s2/L]	0.33	—	[17, 41]
Cve [L / cm H2O]	0.005	estimated from adult	[14]
Rve [cm H2O s/L]	20	estimated from adult	[14]

See Table 1 (Glossary) for variable definitions.

Table 3

Parameters varying with chest wall compliance and simulation conditions.

Parameter	High Cw		Low Cw		Formula	References
	normal Ru	increased Ru	normal Ru	increased Ru
FRC [ml]	24.9	24.9	28.1	28.1	Pel|FRC + Pcw|FRC = 0	[18, 23, 45]
Amus [cm H2O]	1.85	3.2	2.78	3.8	estimated	—
dw [cm H2O]	0.48	0.48	2.4	2.4	estimated	—
Ru,m [cm H2O s/L]	20	200	20	200	—	[39, 41]
Ku [cm H2O s/L]	60	600	60	600	estimated from adult	[14, 39, 41]

See Table 1 (Glossary) for variable definitions.

Table 4

Aggregate parameters and output states targeted during simulations.

Parameter	Ref. Value	High Cw		Low Cw		Formula	References
		normal Ru	increased Ru	normal Ru	increased Ru
CL [ml/cm H2O]	2.3	2.7	2.1	2.3	2.1	(V_A|_EI-V_A|_EE)(P_el|_EI-P_el|_EE)	[8, 39, 40]
Cw [ml/cm H2O]	8.5	9.9	16.0	2.7	3.3	(V_cw|_EI-V_cw|_EE)(P_cw|_EI-P_cw|_EE)	[8, 39]
Crs [ml/cm H2O]	1.8	2.1	1.9	1.2	1.3	(1/CL + 1/Cw)−1	[39, 46]
Rrs [cm H2O s/L]	40	34 to 41	32 to 223	33 to 36	32 to 223	Ru + Rc + Rs	[47]
V˙_E [ml/min]	360	359.2	358.2	360.0	360.0	—	[23]
VT [ml]	6	5.99	5.97	6.00	6.00	—	[40, 48, 49]
V˙_A [ml/s]	±20	-19.4 to 20.8	-16.4 to 28.9	-20.3 to 20.4	-17.1 to 26.4	—	[40, 48]
PA [cm H2O]	±1 − 2	-0.75 to 0.84	-0.96 to 3.63	-0.69 to 0.69	-0.89 to 3.80	—	[23]
Pel [cm H2O]	1 to 6	0.9 to 3.6	2.7 to 5.9	1.8 to 4.8	2.8 to 6.1	—	[40]
Ppl [cm H2O]	-3 to -6	-0.6 to -3.8	0 to -6.2	-1.5 to -4.9	0 to -6.4	—	[23]

EI: end-inspiratory. EE: end-expiratory. See Table 1 (Glossary) for variable definitions.

Computational procedures

All simulations proceeded with an initial respiration rate of 60 breaths/min (f = 1), initial minute ventilation , and initial tidal volume V = 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 regardless of chest or lung compliance; and two, this allowed for comparison of simulation results originating from similar starting points. Distinct values for A were prescribed for each simulation to achieve the initial , 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. CPAP is simulated with an increase of P from 0 to 5 triggered at F = 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 . Each simulation was performed three times: constant f; variable f by breath according to where V 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.

Table 5

Initial conditions.

Initial Condition	High Cw	Low Cw	Formula	References
V˙(0)	0	0	—	—
Vc(0)	0.0001	0.0001	estimated from adult	[12]
Pel(0)	0.954	2.015	Pel|FRC	[18]
Pve(0)	0	0	—	—

Results

Parameterized static compliance curves for V(P) and V(P) are shown in Fig 2 for high C (left) and low C (right). The hysteretic tidal breathing loops are superimposed on the curve V(P) for normal R in black and increased R in grey. Hysteresis is caused in the model by the viscoelastic parameters C and R, which were tuned to maintain appropriately valued lung volume outputs.

Fig 3 shows the impact of high vs low C and normal vs. high R on the five states P, P, P, V, and . Increased R increases P almost threefold, but C has very little impact. However, decreased C increases P significantly, effectively raising it higher on the lung PV curve. Increasing R even higher increases P but there is no difference with respect to C. The opposite appears to occur with P dynamics, in that decreasing C makes P more negative (“increasing” the magnitude of the pressure) and increasing R strengthens that effect. Low C and subsequently high R increase V, mimicking the effect for P. High R shifts 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 H2O, -30 to 30 ml/s, and 8 ml, respectively.

Fig 3

Simulated periodic steady-state tracings of five breaths.

Depicted are alveolar volume, airflow, alveolar pressure, dynamic elastic lung recoil, and pleural pressure, under high and low C conditions, with normal vs. high R.

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 compared to the same with low C. 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 and increased R (S8).

Table 6

Simulations and time to failure (TTF, in hours), defined as 90% volume loss.

Intervention	Cw condition	Variable γ	Simulation	TTF, hours
				constant f	variable f	variable f + AE
None	Low	No	1	2.49	2.53	2.51
		Yes	2	2.25	2.28	2.26
	High	No	3	0.30	0.32	0.30
		Yes	4	0.27	0.29	0.27
Increased Ru	Low	No	5	24.7	24.7	24.5
		Yes	6	22.2	22.2	21.9
	High	No	7	18.5	18.4	17.3
		Yes	8	16.6	16.5	15.0
CPAP, 10% loss	Low	No	9	2.56	2.61	2.60
		Yes	10	2.32	2.36	2.34
	Yes	No	11	0.83	0.79	0.76
		Yes	12	0.83	0.80	0.77
CPAP, 5% loss	High	No	13	2.94	2.50	2.46
1-1 4-7 CPAP, 3% loss			14	8.57	6.89	6.83

Increased R: A 10-fold increase in R was applied during expiration. CPAP: Simulated administration of P = 5 occurred when recruited fraction was down 10%, then again at 5% and 3% with constant γ. AE: A single 20 second apneic event occurred at the 2 minute mark of the simulation.

The breath-to-breath change in EELV and V under high and low C conditions with no interventions are given in Fig 4 (Simulations 1 and 3). The high C simulation reaches accelerated loss of volume and eventually failure at 0.3 hours, much more quickly than the low C 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 4

Breath-to-breath volumes.

End-expiratory lung volume (left y-axis) and tidal volume (right y-axis) under high and low C conditions, no interventions.

Fig 5 shows changes in dynamic lung compliance and tidal volume with high and low C without changes in γ, then adding CPAP to the high C 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 to 5 cm H2O when F < 0.9, which happened when the lung volumes were already decreasing quickly towards failure. However, CPAP triggering at F < 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 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 on dynamics induces the same lung volume loss without other mitigating actions.

Fig 5

Breath-to-breath dynamic lung compliance and tidal volume.

Depicted are high and low C conditions, with simulated CPAP triggered in the high C condition when recruited fraction dropped 10%, 5%, and 3%.

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 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 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 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-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.

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)

Click here for additional data file.

Minimal representative code.

Attached MATLAB code runs Simulation 3 (high C, normal R) with constant frequency for 6 periods.

(PDF)

Click here for additional data file.

Table 1

Glossary.

Parameter/State	Physiologic description
TLC [ml]	Total lung capacity
RV [ml]	Residual volume
FRC [ml]	Functional residual capacity
VC[ml]	Vital capacity
RR [br/min]	Respiratory rate
f [br/s]	Respiratory frequency
T [s]	Duration of respiratory cycle
VT [ml]	Tidal volume
V˙_E [ml/min]	Minute ventilation
V˙_A [ml/s]	Airflow
Amus [cm H2O]	Muscle pressure amplitude
Ptm [cm H2O]	Transmural pressure
PA [cm H2O]	Alveolar pressure
Pel [cm H2O]	Lung elastic recoil (transpulmonary pressure)
Pve [cm H2O]	Viscoelastic component of pressure
Pl,dyn [cm H2O]	Dynamic pulmonary pressure
Ppl [cm H2O]	Pleural pressure
Pcw [cm H2O]	Chest wall elastic recoil
Pmus [cm H2O]	Respiratory muscle pressure
CA [ml/cm H2O]	Lung compliance
Cw [ml/cm H2O]	Chest wall compliance
Crs [ml/cm H2O]	Respiratory system compliance
Rrs [cm H2O s/L]	Respiratory system resistance
ν	Fraction of VC for chest wall relaxation volume
V0 [ml]	Chest wall relaxation volume
β	Baseline fraction of lung recruited at Pel = 0
γ	Maximum recruitable function of lung
α	Lower asymptote, fraction recruitment
k [1/cm H2O]	Characterizes slope, aggregate lung elasticity
cF [cm H2O]	Pressure at maximum lung recruitment
dF [cm H2O]	Characterizes slope at maximum lung recruitment
aw [ml]	Lower asymptote, chest wall compliance
bw [ml]	Characterizes slope, Pcw → ∞
cw [cm H2O]	Transition point, chest wall compliance
dw [cm H2O]	Characterizes slope, Pcw → ∞
ac [ml]	Lower asymptote, collapsible airway
bc [ml]	Upper asymptote, collapsible airway
cc [cm H2O]	Pressure at peak collapsible airway compliance
dc [cm H2O]	Characterizes slope, peak coll. airway compliance
Kc	Collapsible airway resistance coefficient
Vc,max [ml]	Peak collapsible airway volume
Rs,m [cm H2O s/L]	Minimum small airway resistance
Rs,d [cm H2O s/L]	Change in small airway resistance
Ks	Small airway resistance low pressure coefficient
Iu [cm H2O s2/L]	Upper airway inertance
Ru,m [cm H2O s/L]	Laminar value, upper airway resistance
Ku [cm H2O s/L]	Turbulent coefficient, upper airway resistance
Cve [L / cm H2O]	Lung viscoelastic compliance
Rve [cm H2O s/L]	Lung viscoelastic resistance