## Abstract

A computational model for maximal expiratory flow in constricted lungs is presented. The model was constructed by combining a previous computational model for maximal expiratory flow in normal lungs and a previous mathematical model for smooth muscle dynamics. Maximal expiratory flow-volume curves were computed for different levels of smooth muscle activation. The computed maximal expiratory flow-volume curves agree with data in the literature on flow in constricted nonasthmatic subjects. In the model, muscle force during expiration depends on the balance between the decrease in force that accompanies muscle shortening and the recovery of force that occurs during the time course of expiration, and the computed increase in residual volume (RV) depends on the magnitude of force recovery. The model was also used to calculate RV for a vital capacity maneuver with a slow rate of expiration, and RV was found to be further increased for this maneuver. We propose that the measurement of RV for a vital capacity maneuver with a slow rate of expiration would provide a more sensitive test of smooth muscle activation than the measurement of maximal expiratory flow.

- mathematical model
- mechanics
- constricted lungs
- flow-volume curve

some years ago, lambert et al. presented a computational model for maximal expiratory flow from the lung (11). In this model, the areas of the trachea and 16 generations of a symmetric bronchial tree were specified as functions of transmural pressure. For a given lung volume, flow was prescribed, and the convective and dissipative pressure losses were integrated along the airways, from the entrance to the 16th generation to the end of the trachea, to obtain the pressure distribution within the bronchial tree. The calculation was repeated for increasing values of flow until flow speed equaled wave speed at some point in the bronchial tree or the pressure in the flow became more negative than −100 cmH_{2}O. The flow at which one of these conditions was met was identified as maximal flow at the given lung volume. The calculation was repeated for values of lung volume that spanned the vital capacity (VC), and a maximal expiratory flow-volume (MEFV) curve was obtained. The calculated flow-volume curves for air and gases containing He and SF_{6} agreed well with the average curves measured in normal subjects, and the details of the calculated pressure distribution agreed well with the somewhat limited data on the relative contributions of convective and dissipative losses to the total pressure drop during maximal flow and data on the location of the flow-limiting site as a function of lung volume.

More recently, Anafi and Wilson (1) have presented a mathematical model for the dynamic length-tension behavior of airway smooth muscle. This model consists of a set of differential equations relating muscle force, length, and the time derivatives of force and length. Solutions to these equations for a given oscillatory length history match the observed force-length behavior of smooth muscle for different amplitudes of length excursions and a wide range of oscillatory frequencies (3, 20). Anafi and Wilson also incorporated this model into the model of Gunst et al. (5) for the mechanics of an intact airway containing smooth muscle and subjected to a transmural pressure that includes the forces exerted by the attachments to the surrounding parenchyma, the interdependence force. They calculated the time course of muscle shortening and airway narrowing that would occur after a deep breath, and this agreed well with the observed time course of the increase in airway resistance after a deep breath (18).

Here, we describe a computational model for maximal flow in constricted lungs. The model was generated by introducing the Anafi model for constricted airways into the Lambert model for expiratory flow. Flow-volume curves were calculated for different levels of muscle activation. The computed flow-volume curve for maximal muscle activation agrees well with data for flow in maximally constricted normal subjects, and the model provides insight into the mechanics of flow in constricted lungs.

## MODEL

All of the elements of the model have been described in detail in earlier papers, and these descriptions will be reviewed briefly here.

#### Airway geometry and mechanics.

As in the original Lambert model (11), the conducting airways are represented by the trachea (*z* = 0) plus 16 generations of a symmetric dichotomously branching network (*z* = 1–16) extending peripherally from the trachea. Thus each generation consists of 2^{z} identical airways in parallel.

The model for the mechanics of the airways is a combination of the models proposed by Lambert et al. (11) and Anafi and Wilson (1). A schematic diagram of an airway cross section is shown in Fig. 1. The radius and area of the lumen of the airway are denoted *r* and *A*, respectively. The airway wall consists of a layer of incompressible tissue surrounded by a thin band of smooth muscle with radius *r*_{m}. The airway is embedded in the surrounding parenchyma.

Pressure in the lumen of the airway is denoted P. This pressure is balanced by forces exerted by the tissue components of the airway and the surrounding parenchyma. The first of these is the elastic recoil pressure generated by distortion of the connective tissue in the airway wall (Pel). The second is the interdependence stress in the parenchyma that is the result of distortion of the surrounding parenchyma (Pint). The third is the compressive pressure that results from active force in the layer of smooth muscle (Pm). Thus the equilibrium equation for the airway wall is the following. (1) The terms on the right side of *Eq. 1* are functions of *A*, and these functions are described in the following paragraphs.

The relation between Pel and *A* is the same as that used by Lambert et al. (11). (2) The parameters in *Eq. 2*, maximal area (*A*_{max}), normalized area when Pel = 0 (*α*_{0}), P_{1}, P_{2}, *n*_{1}, and *n*_{2}, have different values for each generation in the bronchial tree, and the values that are used here are the same as those listed in Ref. 8.

For the intraparenchymal airways (*z* = 3–16), the relation between Pint and *A* is the one proposed by Lai-Fook (7) and subsequently used in models of airway mechanics (1, 2, 5). (3) In this equation, Pl denotes lung recoil pressure, and *x* describes the mismatch between the outer radius of the airway and the radius of the hole in the undeformed parenchyma. It is assumed that, for the conditions, *A* = *A*_{max}, *r*_{m} equal to its corresponding maximal value, denoted *r*_{0}, and lung volume equal to total lung capacity (TLC), airway radius matches the radius of the hole in the undeformed parenchyma. At other lung volumes, the radius of the hole in the undeformed parenchyma is assumed to equal *r*_{0}*v*^{1/3}, where *v* denotes the ratio of lung volume to volume at TLC. The variable *x* in *Eq. 3* is defined as the difference between the radius of the airway and the radius of the undeformed hole, nondimensionalized by the radius of the undeformed hole, and *x* is therefore given by the following equation. (4) To complete this description of the relation between Pint and *A*, the relation between *r*_{m} and *A* is required. The area of the tissue inside the muscle is assumed to be 16% of the area inside the muscle at *r*_{m} = *r*_{0}. This value for the fractional area occupied by tissue is consistent with the data of James et al. (6). The tissue is assumed to be incompressible, and, therefore, *A* and *r*_{m} are related by the following equation. (5) Finally, the relation between Pm and *A* is the same as that proposed by Gunst et al. (5). The active force in the muscle layer per unit distance in the axial direction along the airway (transverse to the plane of Fig. 1) is denoted F. That is, F equals the active stress in the muscle times the thickness of the layer. The compressive pressure exerted by the muscle layer is given by the ratio of F to *r*_{m}. (6) The value of F depends on the level of muscle activation, muscle length, and length history. The muscle layer is assumed to have adapted to generate maximal force (F_{0}) when *r*_{m} = *r*_{0}. Data for the relation between isometric force (Fiso) and muscle length (4) are approximated by the following equation. (7) During dynamic muscle shortening, F is smaller than Fiso by the factor *f*. That is, (8) Combining *Eqs. 6*–*8* yields the following expression for Pm. (9) The factor F_{0}/*r*_{0} in this equation describes the amount of smooth muscle in the airway and the degree of activation of the muscle. Based on the data of Gunst and Stropp (4), for maximal muscle activation, this factor is assigned the value of 40 cmH_{2}O. Finally, the factor *f* in *Eq. 9* depends on muscle history. For muscle shortening, the factor *f* is determined as the solution to the following empirical differential equation, where *k* and *a* are constants and the superscripted dot denotes the time derivative. (10) The values of *k* and *a* that were obtained by fitting solutions of this equation to the data of Shen et al. (20) are *k* = 28 and *a* = 0.12 s^{−1}. Thus the first term on the right side of *Eq. 10* describes a very sharp decrease in *f* with decreasing *r*_{m}, and the second term describes a relaxation of *f* toward the value 1 with a time constant of ∼8 s.

The relation between *A* and P for airways with active smooth muscle is, therefore, given by the solution to *Eq. 1* for different values of P, using *Eqs. 2*, *3*, and *9*, together with the auxiliary *Eqs. 4*, *5*, and *10* to evaluate the terms on the right side of *Eq. 1*. The relation between *A* and P depends on the history of the maneuver.

An example of a particular *A*-P relation is shown in Fig. 2. This example shows the airway-area curve for the airways in *generation 6* at a value of Pl = 20 cmH_{2}O. In Fig. 2, the value of *A* as a function of Pel is shown by the dashed line. This is the airway-area curve that was used in the original computational model (8). The curve shown by the dotted dashed line is the curve *A* vs. Pel + Pint. This is the airway-area curve with interdependence forces included. Interdependence has little effect near *A*/*A*_{max} = 1, but it has a significant effect for *A*/*A*_{max} small. The magnitude of Pint depends strongly on Pl and hence on lung volume; at lower lung volumes, Pint is considerably smaller. Finally, the curve for P = Pel + Pint + Pm for F_{0}/*r*_{0} = 40 cmH_{2}O is shown by the solid line. Because muscle force depends on the history of the maneuver, the history for this example must be specified. The example shown depicts the *A*-P curve at the moment that maximal flow begins, and the previous history is the inhalation that precedes forced expiration. During that inhalation, the lung and airways expand, and the smooth muscle lengthens. Data show that, during either the lengthening phase of cyclic length changes (20) or lengthening from an initial isometric state (3), force rises slightly above the isometric value. We, therefore, assume that, at the end of inspiration, *f* = 1, and airway area is given by the solution to *Eq. 1* with *f* = 1. For *generation 6* with P = 20 cmH_{2}O and *f* = 1, *A*/*A*_{max} = 0.32 and *r*_{m}/*r*_{0} = 0.66. Thus, for maximal smooth muscle force, airway area is significantly reduced. As flow begins, P drops below 20 cmH_{2}O, the airways narrow, and muscle shortens rapidly. For rapid muscle shortening, the first term on the right side of *Eq. 10* dominates, and the solution to *Eq. 10* is *f* = exp[*k*(*r*_{m}/*r*_{0} − 0.66)]. The curve shown in Fig. 2 is the solution to *Eq. 1* for this relation between *f* and *r*_{m}/*r*_{0}. Thus this curve describes the locus of possible values of *A* and P for the initial flow. Because the force drops rapidly with decreasing muscle length, the curve for Pm ≠ 0 merges with the curve for Pm = 0 as airway area decreases.

### Flow and pressure.

The pressure gradient in the airways is described by the same equation used in the modeling of maximal flow in human lungs (8). (11) In this equation, *s* denotes distance along the airway in the flow direction, Re denotes the Reynolds number of the flow, μ and ρ denote gas viscosity and density, respectively, and V̇ denotes volume flow rate in an airway. The first two terms on the right side of *Eq. 11* describe dissipative pressure losses in the airways, and the numerical coefficients in these terms were obtained from Reynolds's data on pressure losses in flow through a cast of a human bronchial tree (19). The first of these is proportional to flow and gas viscosity, and this term is like the term that describes the pressure loss in laminar flow in a cylindrical tube. The second is like the term that describes the turbulent pressure loss in a smooth cylindrical tube. The third term on the right side of *Eq. 11* describes the convective pressure decrease in flow through a tube with decreasing cross-sectional area.

### Computational method.

The computation of maximal expiratory flow was begun by calculating flow for Pl = 20 cmH_{2}O. The *A*-P curves for this initial flow were calculated as described above. A small value of V̇ was chosen, and the pressure distribution in the bronchial tree was computed by using the same method that was used before (11). First, the Bernoulli drop from the stationary alveolar gas to the flowing gas in the entrance to *generation 16* was calculated. Then *Eq. 11* was integrated along the airway. Then the Bernoulli drop from the end of *generation 16* to the entrance to *generation 15* and the pressure drop through *generation 15* were calculated, and so on through the bronchial tree to the exit of the trachea. In all steps, the dependence of airway area on pressure was taken into account.

The value of V̇ was successively increased, and the calculation was repeated until a value of V̇ was reached for which flow speed equaled wave speed at some point in the bronchial tree, or P at the end of the trachea reached −100 cmH_{2}O. The flow for which one of these two conditions was met was taken as maximal flow at the given lung volume. Flow increments were reduced as maximal flow was approached.

For a small time step, the changes in lung volume and Pl were calculated from the value of maximal flow at Pl = 20 cmH_{2}O, and the change in *f* was computed from the second term on the right side of *Eq. 10*. New *A-*P curves were calculated, and a value of maximal flow for this volume was obtained. These calculations were repeated for successive time steps until a volume was reached for which maximal flow was essentially zero. In this process, the time step was adjusted to be appropriate to the rate of muscle shortening. The initial time step was 5 ms, followed by 20 ms, and 100 ms late in expiration.

## RESULTS

The MEFV curves calculated for different values of F_{0}/*r*_{0} are shown in Fig. 3. For values of F_{0}/*r*_{0} < 20 cmH_{2}O, smooth muscle activation has little effect on flow. For F_{0}/*r*_{0} > 20 cmH_{2}O, maximal flow decreases, and residual volume (RV) increases with increasing F_{0}/*r*_{0}. For maximal smooth muscle activation, maximal flow is less than one-half that in the control state at all lung volumes, and RV is increased by 16% of VC.

Isovolume pressure-flow curves for F_{0}/*r*_{0} = 40 cmH_{2}O, and lung volumes of 75, 50, and 25% VC are shown in Fig. 4. The shapes of these curves are similar to the shapes of the isovolume pressure-flow curves for the control state. The limitation on flow due to turbulent dissipative losses and the wave speed limit on flow depend on essentially the same parameters (9). Therefore, these two mechanisms for flow limitation are inextricably mixed, and, at higher lung volumes, flow limitation is the result of the combination of turbulent and convective pressure losses. At 25% VC, flow is limited by viscous losses, and, at the pressures at which limiting flow is reached, flow speed is significantly lower than wave speed throughout the bronchial tree.

## DISCUSSION

In the modeling described here, a well-validated model for flow limitation in normal lungs and a well-validated model for the dynamic properties of smooth muscle in constricted airways have been combined to obtain a model for maximal flow in constricted lungs. It seems reasonable, therefore, to suppose that the resulting model reasonably represents the constricted state of an otherwise normal lung. However, some questions about the validity of the model remain. The parameters of the model for muscle dynamics were chosen to fit data on the behavior of canine trachealis muscle, and the properties of normal or asthmatic human airway smooth muscle may be different. Furthermore, parallel inhomogeneities are not included in the model, whereas airway constriction may well be nonuniform. Certainly, during quiet breathing, the ventilation distribution in constricted lungs is extremely nonuniform, and it appears that a fraction of the airways are nearly closed (2, 14). However, an inspiration to TLC opens all airways, and the airways would be expected to contract more uniformly during the subsequent fast expiration than during quiet breathing. Nonetheless, due to intrinsic inhomogeneities in airway properties and inhomogeneous muscle activation, inhomogeneities in airway constriction would be expected, and, in fact, regional emptying has been observed to be more inhomogeneous in constricted than in normal dog lungs (12). Therefore, the homogeneous model must be taken to be a description of the average behavior of the constricted lung.

The major results of the modeling are the computed MEFV curves shown in Fig. 3. The relationship between the shapes of the curves for the constricted and control states shown in this figure is like those in examples reported in the literature (e.g., Refs. 10, 16). The most extensive study of maximal expiratory flow in constricted normal lungs is that of Moore et al. (15). These authors measured the dose-response relationship of MEFV curves to inhaled methacholine in 73 normal subjects of both genders and a wide range of ages. Although not all of the subjects reached a plateau in response to increasing doses of methacholine, many did, and one might expect that muscle activation was near maximal at the maximal response. Moore et al. reported the fractional decrease in flow at 30 and 50% VC. The average maximal response of the subjects was a 53% reduction in flow at 50% VC and a 66% reduction at 30% VC. Points with these percentage decreases in flow are shown in Fig. 3. These points lie close to the model prediction for maximal muscle activation. Also, Moore et al. report that the average maximal decrease in forced expiratory volume in 1 s (FEV_{1}) in their subjects was 24%, and, for maximal muscle activation in the model, the decrease in FEV_{1} is essentially the same, 26%. Thus the quantitative agreement between the predicted flows and the data is good. We would like to emphasize the fact that the values of the parameters in the current model were all set by previous matching to independent anatomic and physiological data, and no parameter values were adjusted to match the data for flow in constricted lungs.

A computational model constitutes a quantitative expression of our understanding, and that is its primary use. To the extent that predictions of a model can be compared with data, it provides a quantitative test of our understanding. Also, with a model, quantities that are inaccessible to measurement can be evaluated, and the results of impossible experiments can be calculated, and these exercises can extend our understanding. In the following paragraphs, the information that the model provides about the mechanisms of flow limitation in constricted lungs is discussed.

One notable feature of the curves shown in Fig. 3 is the insensitivity of the MEFV curve to lower levels of smooth muscle activation. For F_{0}/*r*_{0} = 20 cmH_{2}O, the MEFV curve and FEV_{1} are nearly the same as those for the control state. For this value of F_{0}/*r*_{0}, *A*/*A*_{max} for the peripheral airways is 0.7. As a result, the pressure losses in the periphery, which are proportional to *A*^{−2}, are twice those in the control state. However, the peripheral losses in the control state are small, and, with an increase by a factor of 2, they remain small. More importantly, at lower levels of activation, the *A-*P curve for the constricted airway merges with the *A*-P curve for the unconstricted airway at positive pressures, rather than at negative pressures, as shown in Fig. 2. As a result, for negative values of P, the areas of the constricted airways are nearly the same as the areas of unconstricted airways. Thus the *A*-P curve in the region where flow limitation occurs is nearly the same as in the control case. Because FEV_{1} is insensitive to lower levels of muscle activation, one might expect the dose-response curve of FEV_{1} to methacholine (15) to be narrower than the dose-response curve of smooth muscle force in vitro (13), and this is the case.

At higher values of F_{0}/*r*_{0}, this drop in muscle force with decreasing muscle length also diminishes the effect of muscle activation on the MEFV curve. To illustrate this point, an MEFV curve for F_{0}/*r*_{0} = 40 cmH_{2}O, and no effect of length on force (*k* = 0 in *Eq. 11*) was calculated, and this curve is shown in Fig. 5. In this case, muscle force remains equal to the isometric force for maximal muscle activation as the muscle shortens. Flow is drastically reduced at higher lung volumes and drops sharply with decreasing volume. The airways close at a lung volume for which Pl = 10.4 cmH_{2}O. This is consistent with the observation that, under static conditions, maximally activated smooth muscle closes the airways for values of Pl is less than ∼10 cmH_{2}O (5).

The sharp drop in muscle force with decreasing lung volume explains the difference between maximal flows initiated at lung volumes lower than TLC and maximal flows for the full forced VC maneuver (17). For the full maneuver, *f* and muscle force decrease significantly during the volume decrease from TLC to the volume at which the partial maneuver begins. If it is assumed that *f* = 1 at the initial volume for the partial maneuver, *f* and muscle force are higher than for the full maneuver at that volume and remain higher for all succeeding lung volumes. As a result, flows are lower for the partial than for the full maneuver.

The second term on the right side of *Eq. 10*, the equation that governs smooth muscle dynamics, also plays a significant role in determining the MEFV curves shown in Fig. 3. This term describes the recovery of force over time that occurs when muscle force is depressed below the isometric force. To illustrate the effect of this term, an MEFV curve for F_{0}/*r*_{0} = 40 cmH_{2}O and no force recovery (*a* = 0 in *Eq. 10*) was calculated, and this curve is shown in Fig. 5. It can be seen from Fig. 5 that this term affects the MEFV curve at later times in the expiration, and it is the term that determines the magnitude of the increase in RV in the constricted lung.

It can be seen from Fig. 3 that, at F_{0}/*r*_{0} = 40 cmH_{2}O, flow is decreasing rapidly with increasing F_{0}/*r*_{0}. Therefore, muscle hypertrophy would be expected to have a big effect on flow. To model muscle hypertrophy, the MEFV curve for a 25% increase in F_{0}/*r*_{0} from 40 to 50 cmH_{2}O was calculated, and the result is shown in Fig. 6. The increase in F_{0}/*r*_{0} caused a nearly parallel shift of the MEFV curve, and FEV_{1} decreased to 46% of the control value. The effect of inflammation was modeled by increasing the amount of tissue in the wall. Maximal flow was calculated for F_{0}/*r*_{0} = 40 cmH_{2}O and a 25% increase in the tissue volume, from 16% of the area inside the muscle ring to 20%. The result is shown in Fig. 6. The change in wall area also caused a nearly parallel shift in the MEFV curve, but the magnitude of the shift was less than one-half that for the model of muscle hypertrophy. In asthma, muscle hypertrophy and inflammation would be expected to be nonuniform, whereas, in the modeling, the changes in airway properties were applied uniformly. Therefore, these results only indicate the potential magnitude of the effects of muscle hypertrophy and inflammation.

Because RV depends on force recovery and force recovery depends on time, we reasoned that the increase in RV would be greater if expiration were prolonged. We, therefore, calculated the value of RV for a slow expiration. Flow was set at a value of 200 ml/s, and this flow was continued as lung volume decreased until flow limitation reduced flow below 200 ml/s. The resulting curve for F_{0}/*r*_{0} = 40 cmH_{2}O is shown in Fig. 5. The value of RV is higher than the value of RV for the forced VC maneuver. A plot of the increase in RV for this maneuver, expressed as percentage of VC, vs. F_{0}/*r*_{0} is shown in Fig. 7. For smaller values of F_{0}/*r*_{0}, the change in RV with F_{0}/*r*_{0} for slow expiration is considerably greater than the change in FEV_{1} for the forced VC maneuver. That is, RV for the slow-expiration maneuver is more sensitive to F_{0}/*r*_{0} than FEV_{1}, and we suggest the measurement of RV for slow expiration as a more sensitive test of the responsiveness of smooth muscle to agonists.

## 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 © 2005 the American Physiological Society