## Abstract

Ventilator-induced lung injury has been proposed as being caused by overdistention and closure and reopening of small airways and alveoli. Here we investigate the possibility that heterogeneous constriction increases airflow-related shear stress to a dangerously high level that may be sufficient to cause injury to the epithelial cells during mechanical ventilation. We employed an anatomically consistent model of the respiratory system, based on Horsfield morphometric data, and solved for the time evolution of pressure and flow along the airway tree during mechanical ventilation. We simulated constant-flow ventilation with passive expiration in two different conditions: baseline and highly heterogeneous constriction. The constriction was applied with two strategies: establishing a simple diameter reduction or adding also a length shortening. The shear stress distribution on airway walls was analyzed for airways ranging from the trachea to the acini. Our results indicate that *1*) heterogeneous constriction can amplify the maximal values of shear stress up to 50-fold, with peak values higher than 0.6 cmH_{2}O; *2*) the highest shear stress is found in pathways constricted by 60–80%; *3*) simultaneous diameter reduction and shortening amplifies the shear stresses by three- to fourfold, with shear stresses reaching 2 cmH_{2}O; and *4*) there is a range of airways (diameters from 0.6 to 0.3 mm at baseline) that appear to be at risk of very high stresses. We conclude that elevated airflow-related shear stress on the epithelial cell layer can occur during heterogeneous constriction and conjecture that this may constitute a mechanism contributing to ventilator-induced lung injury.

- lung mechanics
- ventilator-induced lung injury
- simulation model

mechanical ventilation (MV) is often required to manage patients with acute respiratory failure (ARF) of different origins, including those associated with abnormalities in respiratory mechanics (31). However, there is considerable experimental evidence that MV itself can also cause or exacerbate lung injury, resulting in so-called ventilator-induced lung injury (VILI) (9, 23, 28). The mechanisms associated with VILI have not yet been fully elucidated. Two mechanisms have received the most attention. One is related to high airway pressures during MV, termed barotrauma, causing excessive overdistension of alveolar structures (26). More recently, it has become evident that injury results more from overstretching the alveolar tissues than from high pressures, and hence it is now termed volutrauma (10). Another potential mechanism is the repeated recruitment and derecruitment of atelectatic lung units (28). In principle, opening a collapsed airway would require very high-pressure drops across the fluid plug blocking the airway (13). During the process of peeling the airway open, high shear stress (τ) and pressure gradients likely cause epithelial damage (3). This form of VILI, termed atelectrauma (28), is, therefore, associated with ventilation at low end-expiratory lung volumes in injured lungs (19) as well as in normal lungs (4). Several recent studies (8, 28, 32) have postulated that VILI also promotes inflammation leading to multiple organ failure (biotrauma). The upregulation of an inflammatory response by mechanical stimuli (32) has prompted several studies to investigate the microscopic mechanical factors that can lead to altered lung function during MV (2, 16, 33).

This study derives a potential new mechanism associated with VILI. Specifically, we investigate the possibility that heterogeneous lung constriction can increase airflow-related τ on the epithelial layer to unexpectedly high levels. Previous modeling studies (14, 21) have shown that the heterogeneity of airway constriction is a crucial determinant of the change in lung resistance and elastance and that inhomogeneous constriction leads to an extremely uneven flow distribution down to the acinar level. Thus mechanical τ, due to high airflow in small airways during heterogeneous bronchoconstriction, may be excessively elevated, perhaps to levels that can damage the epithelial cells during MV. To establish the feasibility of our hypothesis, we advanced a morphometric model of the human lung during heterogeneous bronchoconstriction in MV and predicted τ distributions among the different airway orders.

## METHODS

*Morphometric model.* We adapted a recently developed model for the analysis of the time course of acinar flow and pressure during MV (21). The model is based on the branching network proposed for the human lungs by Horsfield et al. (15), and it is a careful simplification of previously published computational models of lung mechanics in the frequency domain (14, 30). Briefly, the Horsfield airway tree is an asymmetric bifurcating network with 35 orders. Each order has a specific length, diameter, and recursion index, which identifies the degree of asymmetry of that order. *Order 35* corresponds to the trachea. Therefore, the conducting airways are *orders 8–35*, whereas airway *orders 1–7* represent the pulmonary acini (30). In our model, each airway is represented by a combination of laminar resistance and wall compliance (Fig. 1). The resistance (R) of any airway segment *k* is thus described by the classic Poiseuille equation 1 where η is gas viscosity, *l* is airway length, and *r* is airway radius. The airway wall compliance (Cw) is modeled according to Suki et al. (29) and is given by 2 where *h* is airway wall thickness and *Y* is Young's modulus for the specific composition of the airway wall. We neglected the gas compression in the airways, because this is believed to play a minor role at the breathing frequencies (25). To lower the computational burden, we model explicitly only the conducting airways (where convection is the main determinant of ventilation distribution). The pulmonary acini are lumped into equivalent terminal airways leading to equivalent alveoli (21). The alveolar-tissue element is composed of a shunt gas compression compliance for the alveolus and a viscoelastic tissue model (Fig. 1). We set the parameters of the alveolar-tissue compartments equal to the overall pulmonary load, as assessed in mechanically ventilated, anesthetized, paralyzed normal subjects (5). The impedance of the chest wall is defined by a second-order viscoelastic model by using the parameters reported in Ref. 5. We empirically describe the patient-ventilator interactions by modeling the ventilator as a constant-flow generator during the inspiratory phase, whereas the expiratory phase is assumed to be passive; the switch between these phases is time cycled. The endotracheal tube and the expiratory limb of the ventilator are described as Rohrer resistances according to the values reported by D'Angelo et al. (6). According to the human Horsfield orders and recursion indexes (15), the resulting model contains 106,447 airways, ranging from 0.03-mm diameter at functional residual capacity to the trachea.

*Shear stress.* Assuming laminar flow, the τ exerted by the airflow on the airway wall at generation *k* can be calculated, as shown in the appendix (see *Eq. A11*), according to 3 where V̇ is the flow in the airway segment. The V̇ can be calculated from the relation V̇(*k*) = ξP(*k*)/R(*k*), where ξP is the pressure drop across the airway segment. Hence, from *Eqs. 1* and *3*, one has 4 Thus, to simulate the time course of the tangential τ distribution in the tracheobronchial tree, one has to analyze the resistive pressure drops, lengths, and radii in all of the airway generations during MV.

*Simulation protocol.* We simulated constant-flow, volume-controlled MV with physiological tidal volumes (0.7 liter, i.e., 10 ml/kg for a 70-kg subject), and passive expiration. We used a normal respiratory frequency (15 breaths/min) and a prolonged expiration (inspiratory-to-expiratory ratio of 1:3) to mimic the more common pattern of controlled MV used in patients with ARF.

We simulated two different patient conditions: no constriction (baseline) and highly heterogeneous constriction. We modeled heterogeneous constriction in the lung using the same conditions as Gillis and Lutchen (14). Specifically, we applied normally distributed random perturbations to the peripheral airways (diameter <2 mm) with a 20% mean diameter reduction. The standard deviation of the perturbation was 50%. This resulted in an extremely inhomogeneous lung in which there were closed airways randomly occurring throughout the airway tree. We modeled airway closure defining a threshold radius (10% of baseline) so that any further narrowing resulted in a complete closure of the segment. The heterogeneous constriction was applied with and without simultaneous longitudinal airway shortening according to 5 6 where *r*_{c} is the constricted radius, *l*_{c} is the length, μ is the mean percent constriction, CV is the percent coefficient of variation, and *n* is a normally distributed random number with zero mean and unit variance. The 30° angle employed in *Eq. 6* was used according to the mean value reported from Ebina et. al. (11) for airway smooth muscle orientation in a helical pattern around the airways. Thus we either simultaneously applied *Eqs. 5* and *6*, or only *Eq. 5*. Note that, in the former case, the airway shortening is ∼40% less than airway narrowing. In all cases, there were no changes imposed on the tissue properties of the model. The time course of the distribution of τ was calculated and analyzed for airway generations ranging from the trachea to the acini.

## RESULTS

The time course of τ during the respiratory cycle for the healthy lung in a representative generation (airway *order 10* corresponding to a baseline diameter of 0.4 mm) of the Horsfield model is shown in Fig. 2. The corresponding mouth flow and the airway opening pressure are also shown in the same simulated condition. These signals were recorded after several breathing cycles to avoid including any transient effect in the results. From the analysis of the absolute maximal values reached during MV for all of the airways in the tree, we constructed the distribution of peak τ (Fig. 3), where the symbols represent the mean of the peak values for all segments within a generation and the bars represent the corresponding upper and lower limits of the peak values. Here we referred to absolute magnitudes because we are interested in the absolute effect of τ on the epithelial cell layer. We note that, at baseline, the range of maximal τ for a given order is very small. Also, the distribution of maximal τ is distinct for the different airway orders (ranging from 0.0003 to 0.0046 cmH_{2}O). The differences between generations are mainly due to the ratio between the unobstructed length and radius since, from *Eq. 4*, this is the factor modulating τ for a given pressure drop. Minimal τ are found in the trachea (*order 35*) and at the acinar level, where the cross-sectional area is roughly the sum of all of the branching of the terminal airways. Maximal τ peaks are found in *generation 9* and are in the range from 0.0042 to 0.0046 cmH_{2}O. The small intrageneration variations are due to the intrinsic asymmetry of the Horsfield structure that leads to a slightly inhomogeneous flow distribution; this variability can be noted particularly in *generations 12–9*.

Figure 4 shows the distribution of the absolute peak τ vs. Horsfield order in the heterogeneous constriction with the two different constriction strategies. We found a huge increase in the peak τ in the small airways that is augmented 50- to >200-fold, compared with the unconstricted case, depending on whether length shortening occurred simultaneously with airway constriction. It is worth noting that the severe heterogeneity in airway mechanics and random airway closures has caused a major impairment in airflow distribution, leading to high values of τ in the entire airway periphery (where constriction was imposed). Interestingly, the largest values are found in the same generations (*order 10* in the Horsfield structure), regardless of the constriction strategy. The average values of peak τ (symbols in Fig. 4) are increased in all of the small airways from 6- to 10-fold. These quantities, compared with the amplification in the maximal peak values, provide a measure of the extremely skewed distribution of τ. Two examples of the time course of τ during MV in a representative airway order of the Horsfield structure (*generation 10*) are shown in Fig. 5. Comparing Figs. 2 and 5, we note that the severe heterogeneity in a compliant walled airway tree not only produces a huge increase in τ peak values but also markedly changes the time course of the analyzed traces. This is reflected in both the amplitude and timing of the τ minima and maxima. These results are very sensitive to whether airway constriction occurs simultaneously with airway shortening (Fig. 5). In fact, concurrent longitudinal and circumferential airway compression yields a marked increase (3.1-fold in this airway generation) in the magnitude of peak τ from circumferential compression alone.

## DISCUSSION

MV is a life-saving procedure that is required to support patients whose respiratory function has been acutely compromised by postoperative complications, trauma, exacerbation of chronic lung diseases, or acute respiratory distress. However, there are a number of clinical studies and animal models demonstrating that mechanical ventilators can induce or exacerbate lung damage (10, 19, 22, 26, 27, 33, 34). This has been attributed to macroscopic mechanisms such as overdistention of lung units or repeated opening and closing of airway and alveolar duct units. Recently, there has been increased attention to the possible cellular mechanisms that may promote a ventilator-induced lung inflammatory response and eventually a systemic inflammatory response (32) without gross mechanical damage to the lung. Nevertheless, the relationship among mechanical stimuli, lung injury, and cellular inflammatory response is not clear.

The purpose of this work was to present the following new hypothesis: mechanical τ on airway walls due to airflow during heterogeneous bronchoconstriction is significantly amplified. Consequently, this flow-related τ may be a possible signal for mechanotransduction in the epithelium of the airways. We examined this τ amplification by advancing a previously proposed morphometric modeling approach (21) to predict the time course of τ in the airway tree during heterogeneous peripheral airway constriction and constant-flow MV. Before discussing the biological relevance of our findings, we first examine several factors that can influence our estimates of the magnitude of the τ.

*Mechanisms influencing shear stress.* Which are the airways in which highest τ occur? From *Eq. 3*, the time course of the τ is directly proportional to the flow in any given airway segment. However, *Eq. 3* also indicates that different τ values result from the trade off between the flow (which, at any bifurcation, is governed by the downstream impedance) and the actual airway radius. Figure 6 shows the relationship between airflow, airway radius, and peak τ amplification in a representative airway order of the Horsfield structure (*generation 10*). Note that closed segments are not reported in Fig. 6. Also, the data are normalized with respect to the baseline (unconstricted) simulation, so that Fig. 6 represents τ amplification due to heterogeneous constriction. This inhomogeneous constriction results in a wide range of different airflows. Recall from *Eq. 3* that there is a need for a reasonable flow through a small diameter to greatly increase airway τ. In fact, Fig. 6 shows that the highest τ amplifications with respect to baseline (≥50-fold) occur in quite highly constricted airways (≥60–70%) and with a broad range of peak flows (from -80 to +100% with respect to baseline).

Our model simulates peripheral airway constriction without accounting for parenchymal interdependence. In reality, when an airway constricts, the neighboring airways may be kept, to some extent, open due to parenchymal interdependence (20). For example, if we consider a bifurcation, extreme constriction of one of the daughters may limit the constriction of the other daughter. As a consequence, the majority of the flow from the parent will go through the less constricted airway. Hence, the less constricted airway may see a large flow, which can result in a high τ. Thus parenchymal interactions can potentially further amplify τ in certain airways. In other words, the heterogeneous constriction pattern used in our model is consistent with few occurrences of this type of high τ amplification, but airway-parenchymal interdependence could increase the number of such patterns.

Our simulations indicate that a simultaneous airway constriction and length shortening significantly amplify the τ in the entire periphery of the lung (Figs. 4 and 5). The constriction strategy producing simultaneous length shortening and airway diameter reduction was motivated by the helical arrangement of the smooth muscle within the airway wall. A theoretical study (1) has shown that the orientation of smooth muscle around the airway could have a major impact on airway constriction, influencing both airway radius and length. The τ in our simulations increases significantly with airway shortening, which can be interpreted as an augmented heterogeneity of the airway pressure drops. *Equation 4* suggests that, during our simulations with simultaneous airway narrowing and length shortening, the only term that can vary from segment to segment within a generation is ξP(*k*), i.e., the pressure drop across the airway segment in generation *k*. Because the constricted radius-to-length ratio is constant (see *Eq. 6*), the variability in τ must be due to the variability in ξP(*k*). This means that the pressure drops across certain airway segments must also increase significantly to generate the τ amplification shown in Fig. 6. Interestingly, if we leave out airway shortening and keep just airway length variability in *Eq. 6*, we found that (data not shown) length variability accounts for a large part of the amplification in τ compared with airway narrowing alone.

Another phenomenon not accounted for in our model is that the airway walls are coated with a thin liquid layer, which, in principle, could also influence the local values of the τ in individual airways. However, as detailed in the appendix, the τ at the air-liquid interface is directly transmitted to the airway wall without attenuation. This is true under the assumptions that the fluid is Newtonian and the flow in the liquid layer is approximately plane flow. The assumption of Newtonian fluid may not be valid, especially during injury due to the presence of surfactant and membrane-related lipids and various proteins in the liquid layer. In any case, there is also experimental evidence that the airway surface liquid is quite heterogeneously distributed in the peripheral airways and that, in virtually all airways, both in the circumferential and the longitudinal directions, there are regions where no liquid can be found (36). These regions would be directly exposed to τ without the influence of the surface film.

Additionally, airway walls are not smooth and have ridges and folds; the ratio of the height of the ridges to the diameter will be amplified during constriction (35). The τ in a nonsmooth walled tube are amplified, depending on the ratio of the height of the ridge and the airway diameter (12). In fact, the fluid mechanical calculation of Gaver and Kute (12) suggests a threefold increase in τ in a narrow tube with a ridge-to-tube diameter ratio of 0.1. Because we estimated similar ridge-to-diameter ratios in the guinea pig airway from images published by Yager et al. (36), it is likely that an additional amplification of τ could occur along these ridges within the airways. Thus it is possible that even higher τ amplifications (up to 600-fold) exist in the lung, which are highly localized within certain airways, and peak τ may reach values higher than 2 cmH_{2}O.

Another possible source of τ amplification not accounted for in our model is turbulent flow. *Equation 3* holds only under the conditions of laminar flow. However, at the flow rates used in our simulations, turbulent τ would only exist in the first three to four generations of the model (trachea and main bronchi, where peak Reynold's numbers exceed 2,000) and not in the small airways (that have peak Reynold's numbers of <900, even in heterogeneous constriction). In our model of wall τ, we also neglected the so-called “entrance region” effects, due to airway bifurcation. In the entrance regions, the flow velocity profile is flat, resulting in higher τ at the airway wall, because the velocity gradient near the wall increases (24).

*Significance of shear stress.* It is important to examine whether our calculated τ have any physiological relevance. Our results indicate that heterogeneous constriction markedly amplifies the mean values of τ in the small airways (ranging from 5- to 17-fold). More importantly, however, we found that the amplification of peak values of τ can be much higher than that of the mean τ, reaching absolute τ values of >2 cmH_{2}O. To the best of our knowledge, this is the first time that a value for the τ due to airflow on the airway walls has been estimated. Therefore, it is difficult to predict whether these stresses would be enough to promote inflammation or injury of the epithelial cells. However, it should be emphasized that, during MV, the amplified τ are repeated at a rate of ∼20,000 cycles/day. It is also worth noting that, unlike the endothelial cells, which are designed to tolerate higher blood flow-related τ, epithelial cells may tolerate less τ and hence may be more at risk.

Blood flow-induced τ in the vascular system have been studied and are known to cause endothelial cell remodeling, as well as endothelial mechanotransduction (7, 18). In hemodynamics, τ values ranging from 0.07 to 0.1 cmH_{2}O are considered high; the range of 0.4–1.5 cmH_{2}O is considered supraphysiological. By correcting our τ predictions for the Gaver and Kute effect (12), at least 5% of the airways between orders of 13 and 8 in the Horsfield structure (that is ∼2,550 segments) experience τ peaks >0.24 cmH_{2}O. Therefore, it is possible that these stress levels may play a role also in epithelial cell damage, inflammation, and perhaps even remodeling. Indeed, it has recently been reported that pressure gradients and fluid τ can cause the membrane disruption of epithelial cells in an in vitro model of airway reopening (3). Bilek and colleagues (3) suggested that the major determinant of the number of dead or injured cells is likely to be the magnitude of pressure gradient, especially in conjunction with large τ. The authors measured τ and pressure gradient magnitudes that are much larger than what we have found in our simulations (10- to 20-fold higher). However, they have studied injuries caused by a single reopening event, whereas we believe that, in airflow-related τ, repeated cycles can further amplify the effects of abnormally high stresses.

We believe that the present results may be relevant for ARF conditions in chronic or acute respiratory insufficiency, where highly heterogeneous pathologies exist. Indeed, amplified τ might even be relevant in the triggering of inflammatory or remodeling responses associated with asthma and chronic obstructive pulmonary disease, pathologies also associated with highly heterogeneous constriction. It should be emphasized that, whereas in asthma and chronic obstructive pulmonary disease airway smooth muscle shortening generates both diameter reduction and length shortening, edema and acute respiratory distress syndrome are more governed by surfactant deficiency and alveolar flooding. Therefore, in the latter cases, simultaneous longitudinal and circumferential airway compression are unlikely. Hence, according to our model analysis, we expect a less extreme, but perhaps still significant, amplification of airflow-related τ.

In our simulations, we have only imposed one pattern of inhomogeneous constriction characterized by a relatively high standard deviation of airway narrowing. Obviously, there are an infinite number of possible constriction patterns that can be imposed on the model. However, recent experimental data (17) support the crucial role played by both the mean airway caliber reduction and the heterogeneity of this constriction in asthma, even at baseline. Moreover, these data indicate that asthmatic subjects with severe baseline constriction present frequency dependence of dynamic resistance and elastance consistent with greater heterogeneity than that considered here.

In summary, in the present study, we have employed an anatomically consistent mathematical model of the mechanically ventilated patient and predicted that heterogeneous constriction produces a huge amplification of the τ values on the airway walls compared with those encountered in the unconstricted baseline condition. Conservative estimates predict that local peak τ can reach 2 cmH_{2}O in certain small airways. Because these numbers are considered supraphysiological stresses for endothelial cells, our predictions lead to the hypothesis that physical forces, such as cyclic airflow-related τ due to inhomogeneous ventilation distribution in the small airways, can potentially be an important contributor to the mechanical signals associated with initiating or exacerbating airway inflammation and injury during MV. Nevertheless, experimental studies are needed to assess the biological relevance of the proposed mechanism.

## APPENDIX: SHEAR STRESS DUE TO AIRFLOW IN A CIRCULAR TUBE COVERED BY THIN LIQUID FILM

Let us take a tube of radius *r* that is covered inside by a thin liquid film of thickness *h*. The thickness of the air core is *a*, so that *r* = *a* + *h*, where *h* << *a*. There is a constant pressure gradient across the tube so that air is flowing from left to right. We assume that both the airflow and liquid flow can be considered steady, incompressible, and viscous. We also neglect the effects of gravity and assume circular symmetry. In this case, the Navier-Stokes equations in polar coordinates reduce to a single equation in the axial (*z*) direction A1 where *x* is the radial coordinate, P is the pressure, *v* is the velocity in the *z* direction, and η is the viscosity. The boundary conditions are that *1*) the fluid flow velocity is zero at the wall *v* (*a* + *h*) = 0, and *2*) the τ along the air-liquid interface is continuous. Additionally, for *x* < *a*, the viscosity has a value of air viscosity (η_{a}), whereas, for *a* < *x* < *a* + *h*, the viscosity is equal to that of the liquid (η_{l}). We first calculate the velocity in the air, *v*_{a}. The term in the parentheses of *Eq. A1* can be written as Because the pressure gradient dP/d*z* is constant, we can integrate *Eq. A1* A2 Because the solution has to be finite at *x* = 0, the constant *c*_{1} = 0. Next, we calculate the velocity in the liquid, *v*_{l}. The first and second terms in the parentheses of *Eq. A1* can be approximated as *v*/*h*^{2} and *v*/(*xh*), respectively. Because we assumed that *h* << *a*, the radius of curvature in the liquid layer is large, *v*/*h* >> *v*/*x*, and we neglect the second term in the parentheses of *Eq. A1*. Additionally, we assume that the pressure gradient across the tube is created in the air compartment, and we also drop the dP/d*z* term from the solution for liquid flow. Thus *Eq. A1* is now simplified to A3 The solution of *Eq. A3* is a plane Couette flow A4 The boundary condition at the wall requires that A5 The τ is defined as the product of viscosity and the gradient of strain given by A6 At the air-liquid interface, the τ must be continuous, which implies that A7 The constant *c*_{3} can be obtained from *Eq. A7* and used in *Eq. A5* to obtain *c*_{4}. Inserting these constants into *Eq. A4*, we obtain the velocity solution in the liquid layer A8 Using *Eq. A6* and replacing *a* with *r* - *h*, the final expression for the wall τ is A9 This expression suggests that the pressure gradient is directly balanced by the wall τ. It is interesting to note that the wall τ does not explicitly contain the liquid viscosity. However, the wall τ implicitly depends on the velocity and viscosity of the air through *Eq. A2*. Thus the pressure gradient imposed on the tube drives air through the tube, which in turn generates τ at the air-liquid interface. Under the assumptions we used, the liquid flow can be considered Couette-like, and the τ at the air-liquid interface due to airflow is transmitted to the wall without attenuation.

Note that, under steady laminar flow conditions, fluid inertia can be neglected, and the pressure variation along the tube is constant and equal to pressure drop due to the Poiseuille resistance of the pipe times airflow (V̇) A10 Replacing *Eq. A10* into *Eq. A9*, one obtains the equation for τ at any time point in all of the airway segments A11

## Acknowledgments

The work of G. Nucci was supported by Ministero dell'Istruzione, dell'Università e della Ricerca Scientifica of Italy, Rome. The work of B. Suki and K. Lutchen was supported by National Science Foundation Grants BES-0114538 and BES-0076818 and by National Heart, Lung, and Blood Institute Grant HL-62269-02.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2003 the American Physiological Society