Patterns of recruitment and injury in a heterogeneous airway network model

In respiratory distress, lung airways become flooded with liquid and may collapse due to surface-tension forces acting on air–liquid interfaces, inhibiting gas exchange. This paper proposes a mathematical multiscale model for the mechanical ventilation of a network of occluded airways, where air is forced into the network at a fixed tidal volume, allowing investigation of optimal recruitment strategies. The temporal response is derived from mechanistic models of individual airway reopening, incorporating feedback on the airway pressure due to recruitment. The model accounts for stochastic variability in airway diameter and stiffness across and between generations. For weak heterogeneity, the network is completely ventilated via one or more avalanches of recruitment (with airways recruited in quick succession), each characterized by a transient decrease in the airway pressure; avalanches become more erratic for airways that are initially more flooded. However, the time taken for complete ventilation of the network increases significantly as the network becomes more heterogeneous, leading to increased stresses on airway walls. The model predicts that the most peripheral airways are most at risk of ventilation-induced damage. A positive-end-expiratory pressure reduces the total recruitment time but at the cost of larger stresses exerted on airway walls.

7. The network is mechanically ventilated by injecting a finger of air with a prescribed flow rate. The time-dependent airway pressure is continuous across each airway bifurcation.
8. The influence of surfactant transport on the gas-liquid interface is ignored and so the surface tension coefficient is constant.
9. Of the two steady modes of recruitment elucidated for a single airway, known as 'pushing' (where the air finger pushes a long column of fluid ahead of itself like a leaky piston) and 'peeling' (where the airway walls peel apart rapidly) (Halpern et al., 2005), only 'peeling' motion can lead to appreciable airway recruitment. The recruitment speed can be calculated as a function of the airway pressure, derived in a manner identical to Jensen et al. (2002).
10. Recruitment of a particular airway takes place while the airway pressure exceeds a yield pressure, which is a function of the airway elastic properties. Once the airway pressure has fallen below this threshold pressure, recruitment in that airway halts until the airway pressure increases beyond the threshold again.
11. The generation of airways at the distal end of the network is assumed to open into a compliant compartment representing the airways in the respiratory zone and the alveoli. The relationship between the airway pressure and volume of this compartment is based on that derived by Salazar & Knowles (1964).
We now explain how these assumptions are expressed mathematically.

Model for recruitment of a bifurcating airway network
We consider generations 11-16 of a human lung using 'Model A' from Weibel (1963) which prescribes regular dichotomous branching of generations (Assumption 1). Fully three-dimensional simulations of individual airway recruitment have shown that the surface-tension-driven collapse of a fluid-lined tube is strongly asymmetric (Heil, 1999a,b;Hazel & Heil, 2003), where the airway typically adopts a ribbon-like configuration downstream of the tip of the air finger (see Fig. S1a). To mimic a buckled configuration, in this model each airway is modelled as a planar symmetric flexible-walled channel in cross-section (see Fig. S1b), characterised by a transverse 'height' and effective 'width' in the planes orthogonal to the axis of the channel (Assumption 2). Airway k in generation j is denoted by subscript j and superscript (k). The notation for the present model is summarised in Table. S1. Weibel (1963) (see also Weibel, 1991) measured the lengths and diameters of individual airways and total cross-sectional area of each generation using casts of human lungs taken at 75% total lung capacity (TLC). The mean length and diameter of airways in generations 11-16 and their total cross-sectional area at 75% TLC are listed in Table S2. Following Lambert et al. (1982), we consider this total cross-sectional area as the 'maximal' airway area for all the airways at a given generation, denoted as A m,j .

Estimation of the maximal airway cross sectional area
The mean maximal cross-sectional area for an individual airway in generation j, denoted a m,j , is assumed to take the form a m,j = A m,j 2 −j . (1) Hence, we can write the total cross-sectional area of generation j at any time t as To account for natural variability between airways across a generation we sample the maximal cross-sectional area for each individual airway (a  Figure S1: Schematic of recruitment of a single airway: (a) three-dimensional sketch of airway recruitment showing the buckled 'ribbon-like' configuration Hazel & Heil (2003); the crosssection shape is characterised by and effective 'width' z and a transverse 'height' 2h. (b) cross section of a single airway, resembling 2D flexible walled channel models Halpern et al. (2005); (c) the relationship between airway pressure and recruitment speed for steady recruitment of a single airway.
(a m,j ). In particular we define the relative stiffness parameter s m,j /a m,j and sample its value from a normal distribution with unit mean and variability d (Assumption 4) Quantity Symbol Global properties P (t) airway pressure p airway liquid pressure Q prescribed flow rate of air finger γ uniform surface tension µ uniform liquid viscosity P 0 threshold opening pressure of alveoli (TOP) T 0 constant longitudinal tension applied across the network Properties of generation j A j total cross-sectional area of airways A m,j maximal total cross-sectional area of airways α j ratio of A j to A m,j a m,j maximal cross-sectional area of individual airways h m,j mean maximal transverse displacement of cross-section for airways α 0j ratio of equilibrium cross-sectional area to maximal α ′ 0j gradient of pressure-volume curve at zero transmural pressure (equilibrium) T j longitudinal tension of airways Properties of airway k in generation j U  m,j /π) 1/2 . In addition, we prescribe the effective 'width' of each airway across the tree as   Lambert et al. (1982) proposed a relationship between transmural pressure across the airway wall and total airway cross-sectional area (like a 'tube law') for generations 0-16 of a human lung (see also Lambert & Beck, 2004). They modelled the expansion and collapse of each generation of airways based on anatomical data (Weibel, 1963, see Table S2). Each generation is characterised by an equilibrium area (at zero transmural pressure) and a maximal area (arising at infinite transmural pressure), recognising that airways stiffen when strongly inflated. The total collective equilibrium (or undisturbed) cross-sectional area of a particular generation is denoted as α 0,j A m,j (where α 0,j is also listed in Table S2). In particular, this model describes the relationship between α j , the total cross-sectional area of all the airways at generation j normalised on the the maximal cross-sectional area for that generation (α j = A j /A m,j ), and the transmural pressure P . They used two hyperbolae to describe the inflation and collapse of each generation, appropriately matched together at P = 0. For all the airways at generation j

Model for airway wall elasticity
where n 1j and n 2j are constants (tabulated in Table S2) and where α ′ 0j represents dα j /dP as P → 0 ± . The relationship between α j and P is plotted in Fig. S2(a) for generations 11-16. Similarly, the pressure-area curve for an individual airway in generation 14 with maximal area a m,j is shown as a solid line in Fig. S2(b). Also shown are the pressure-area curves for airways with maximal area (1 ± 0.1)a m,j . The corresponding compliance of this airway is the gradient of the P against a j curve evaluated at P = 0. This figure illustrates that an increase in the maximal area of an airway corresponds to a decrease in compliance.
Using this notation, the equilibrium transverse 'height' of the cross-section in an individual airway can be expressed h ensuring (from (4)) that z m,j (no sum over j).  Figure S2: The variation in α j , the total generation cross-sectional area A j normalised on the 'maximal' cross-sectional area A m,j , with P , the airway pressure, for generations 11-16. Based on data taken from Lambert et al. (1982), listed in Table S2.
Representing each airway as a two-dimensional channel, Lambert's tube law for generation j (5) can be re-expressed in the form where h (k) j is the transverse displacement of the cross-section for airway k in generation j (Assumption 5).
In addition, a constant longitudinal wall tension T 0 is applied between the upstream and downstream ends of the network. This tension is assumed to be constant for all the airways across a given generation, denoted T j , and is calculated using a simple force balance around each bifurcation. If η is the symmetric branching angle, this implies T j = 2T j+1 cos(η). Assuming small branching angles (η ≪ 1) and neglecting external tethering forces, T j ≈ 2T j+1 .

Model for airway collapse
We assume that each airway cross-section (Fig. S1) constitutes two (finite-length) flexible sheets confining a region of Newtonian liquid with viscosity µ (Fig. S1). Initially each constituent airway is collapsed and each sheet is uniform, held a constant distance of 2H (k) j apart by the liquid (Assumption 6). For simplicity we place an axis midway between (and parallel to) the undeformed sheets, about which we assume the motion is symmetric. Throughout this work we assume that the H (k) j is 100C% of the equilibrium transverse 'height' of the airway wall (h (k) e,j ), where C is a constant. Thus, 0 < C ≤ 1 measures the degree of initial airway collapse. In the main paper we investigate the interval 0.1 ≤ C ≤ 1. The coordinate x measures distance along the airway.
A finger of air is introduced at the upstream end of the network with a prescribed flow rate Q (Assumption 7); the corresponding airway pressure is denoted p b (t) where t is time. The finger of air inflates the airway walls, forming an interface with the liquid with constant surface tension γ. This surface tension coefficient is assumed constant, so the effects of surfactant transport through the network are ignored (Assumption 8). Neglecting the pressure drop due to the surface tension of the liquid (assumed to be much less than the wall tension, Halpern et al. (2005)) implies that p b = P . Henceforth we denote the airway pressure as P (t).

Steady recruitment
When the airway pressure exceeds a threshold for an individual airway, denoted P * j (k) , steady recruitment proceeds through one of two modes, known as 'pushing' (where the air finger pushes a long column of fluid ahead of itself like a leaky piston) and 'peeling' (where the airway walls peel apart rapidly) (Halpern et al., 2005), as illustrated schematically in Fig. S1(c). Previous work considering the recruitment of a bifurcated airway has shown that for sufficiently high flow rates the system can sustain fast quasi-steady peeling motion in both daughter airways. For lower flow rates (where simultaneous steady peeling of both daughter airways cannot be sustained) the finger tip in one or both of the daughter airways makes the transition to slow pushing motion once the pressure falls below the critical pressure required to sustain steady motion. For simplicity, since peeling speeds greatly exceed pushing speeds, we assume that the finger of air exhibits quasi-steady peeling in any daughter airway where the pressure exceeds this critical threshold pressure (Assumption 9). Once the airway pressure falls below this threshold in a particular airway, we assume that recruitment of that airway stops until the airway pressure rises above this threshold again (Assumption 10).
We follow Jensen et al. (2002) in constructing a scaling approximation for the recruitment speed in peeling motion which can be implemented for each airway in the network (Assumption 9), dividing the flow into three asymptotic regions (Fig. S1b).
In region I we assume that P ≥ 0 and the airway wall is inflated to a transverse 'height' h (k) I,j (dependent on P ) and in local equilibrium (we assume h m,j ). Modifying (7b), we express the normal stress balance across the airway wall as combining Lambert's tube law with a linearised expression for the effect of airway wall tension. The transverse 'height' of the airway wall far upstream of the tip of the air finger, where h (k) j is independent of x, is given by Across region I we balance the elastic restoring force with the airway wall tension, represented by the two terms on the RHS of (9). This balance implies that the scale of variation within region I can be expressed approximately as from which we can estimate the peeling angle (Fig. S1b) Region II contains the tip of the air finger and the flow field is complicated, but Halpern et al. (2005) have shown that it can be assumed passive provided the flux of fluid passing the tip and the corresponding pressure drop are accounted for.
Across region III (of approximate length L (k) III,j ) the airway wall collapses to a transverse 'height' H (k) j downstream. We choose this to be a fixed fraction, C, of the equilibrium transverse displacement of the airway, so H e,j . Throughout this work we assume that C is a constant across the network (Assumption 6). Following Jensen et al. (2002), we estimate the pressure gradient in the long expanded region ahead of the finger tip from lubrication theory as where p is the liquid pressure in region III. The airway wall tension is assumed large enough to dominate the negative transmural pressure induced by bending of the airway wall, implying so a second estimate of the peeling angle θ Equating the two estimates (11,14) of the peeling angle θ (k) j gives an expression describing steady 'peeling' speed for airway k in generation j, relating the airway pressure P to the speed of the recruitment U (k) j : where h (k) I,j depends non-linearly on P . However, P must exceed a threshold in each daughter airway, which we denote P ⋆ jk .

The yield pressure for airway recruitment
A suitable description of the speed of airway recruitment in steady pushing motion (denoted U (k) j ) is required to derive an estimate of the critical pressure for steady recruitment in an individual airway P * j (k) (where the steady pushing and peeling branches intersect, discussed Halpern et al. (2005) and illustrated schematically in Fig. S1(c)). The upstream transverse displacement of the airway wall can again be approximated using (9), so the total flux of liquid lost across the meniscus (according to Bretherton, 1961) takes the form whereŨ (k) j is the pushing speed and α s ≈ 1.337. The total volume of fluid swept up ahead of the advancing air finger tip takes the form For steady pushing these two fluxes balance exactly and thus we derive that The critical pressure required for recruitment a particular airway occurs when steady pushing and peeling intersect, so P ⋆ jk is a solution of calculated numerically using Newton's method implemented in MATLAB. For the simulations shown we only accepted solutions which MATLAB returns as converged within a tolerance of 10 −8 . In these simulations we used an initial pressure of 0.5545 cm H 2 0 and an initial speed of 0 cm s −1 for each airway, but found that the solver was not sensitive to the choice of initial conditions. We performed a variety of numerical experiments starting with a wide range of randomised initial conditions; simulations which MATLAB reported as converged always attained the same value. The critical recruitment pressure for each generation using our benchmark parameters (Table S3 below) is shown in Fig. 1(d) in the main paper.

Model for alveolar compliance
We assume that the airways in the most distal generation considered (generation 16 in this case) open into a compliant bag representing the airways in the respiratory zone and the alveoli (Assumption 11). Following Bates & Irvin (2002), we assume that this sac will expand and contract with changing P , according to the Salazar-Knowles relationship (Salazar & Knowles, 1964), where V A is the maximal volume of each acinar compartment and K A is a compliance parameter Bates & Irvin (2002). This model has zero acinar volume for P = 0. We modify this model to allow the alveoli to accommodate a theoretical opening pressure (TOP) for the alveoli P A such that

Ventilating the airway network
As in previous compartmental models (Halpern et al., 2005), we assume that in each airway being recruited the rate of lengthening of the air finger is equal to the recruitment speed identified in (15), so that d dt L where l is an index over all the acini and the function H l is a Heaviside function, which is set to zero for P < P A and set to one for P ≥ P A . Eqn. (23) can be rearranged to form an ODE for dP/dt in terms of P , U j , L (k) j and the model parameters.

Simulating the model numerically
For M airways being recruited at time t, we have a system of 2M equations ( (15) and (22) for each airway for the M quasi-steady recruitment speeds and M lengths of bubble fingers). These are coupled to the pressure condition (23) forming a system of 2M +1 equations. These governing equations were solved using MATLAB solver ode15s for differential-algebraic equation systems. A typical simulation for a 6 generation network takes approximately 2 minutes on a desktop computer. In simulations we employed an absolute tolerance of 10 −6 and a relative tolerance of 10 −4 . Selected simulations have been validated against equivalents using error tolerances a factor of 10 smaller (everything else held the same) and the timetraces of airway pressure were almost indistinguishable.
To complete the description of the numerical model we specify an initial condition for the airway pressure P (t). In this study we explored two possible initial conditions for the system.
The first choice was to assume that theme that the bubble is peeling steadily in the parent generation, with the airway pressure and peeling speed determined directly from the choice of prescribed volume flux, calculated from (19) for j = 11 using Newton's method (with an absolute tolerance of 10 −8 and a relative tolerance of 10 −8 ).
Alternatively, the second choice was to prescribe a positive end expiratory pressure (PEEP) during recruitment, which maintained that the airway pressure must always exceed a particular value P E ; in simulations of the model including a PEEP we prescribed the initial pressure as P = P E .