## Abstract

We do not yet have a good quantitative understanding of how the force-velocity properties of airway smooth muscle interact with the opposing loads of parenchymal tethering and airway wall stiffness to produce the dynamics of bronchoconstriction. We therefore developed a two-dimensional computational model of a dynamically narrowing airway embedded in uniformly elastic lung parenchyma and compared the predictions of the model to published measurements of airway resistance made in rats and rabbits during the development of bronchoconstriction following a bolus injection of methacholine. The model accurately reproduced the experimental time-courses of airway resistance as a function of both lung inflation pressure and tidal volume. The model also showed that the stiffness of the airway wall is similar in rats and rabbits, and significantly greater than that of the lung parenchyma. Our results indicate that the main features of the dynamical nature of bronchoconstriction in vivo can be understood in terms of the classic Hill force-velocity relationship operating against elastic loads provided by the surrounding lung parenchyma and an airway wall that is stiffer than the parenchyma.

- bronchial responsiveness
- computational model
- force-velocity behavior
- force-length behavior

the conventional approach to assessing airways hyperresponsiveness in animal models is to administer an airway smooth muscle (ASM) agonist to the lungs and measure the resulting change in lung function (46). Bronchoconstriction is a dynamic event, however, so degree of airway narrowing is not the only consequence of bronchial challenge that bears on airway responsiveness. Another key issue is how rapidly the airways narrow from baseline to the point of maximal constriction (6–8, 50). Indeed, in some studies, the rate of ASM shortening has been implicated as a key difference between hyperresponsive animals and their controls (13, 26). Nevertheless, the link between airway responsiveness and both rate and magnitude of ASM shortening is complicated by the fact that shortening velocity exhibits a strong inverse dependence on the magnitude of any opposing mechanical load. One such load is provided by the outward pull of the parenchyma in which the airways are embedded, as evidenced by the observations that airway responsiveness is highly inversely dependent on both tidal volume (Vt) (38) and mean lung volume (7, 12). There is also evidence to suggest that the stiffness of the airway wall itself presents an additional load against which the ASM must shorten (1).

We do not yet have a good quantitative understanding of how the force-velocity properties of ASM interact with the opposing loads of parenchymal tethering and airway wall stiffness to produce the dynamics of bronchoconstriction. This issue is clearly important for our understanding of the nature of airway hyperresponsiveness itself, yet has so far received very little attention at the level of the whole animal from the mathematical/computational modeling community. We therefore developed a 2-dimensional computational model in which we represent the entire airway tree as a single “effective” cylindrical airway embedded in uniformly elastic lung parenchyma. The active hoop stress in the airway wall is linked to both airway radius and transpulmonary pressure (Ptp) using equations from the literature. We calculate the dynamics of airway narrowing by integrating the changes in active wall stress and radius over time using the classic Hill force-velocity relationship for smooth muscle. Our goal was to determine the extent to which this simple model can capture the dynamics of bronchoconstriction, as represented by published measurements of airway resistance in rats and rabbits following challenge with an intravenous bolus of methacholine.

## METHODS

#### Model development.

We model an airway in two dimensions as a circular ring of ASM wrapped around an elastic airway wall embedded in homogeneously elastic lung parenchyma. When the ASM contracts, it narrows the airway against the outward pull of the parenchyma that is attached to the outside of the airway wall. To calculate how ASM activation narrows this airway, we must identify an expression for the balance of forces that determine airway radius. These forces are of two types: *1*) passive forces due to the elastic recoil of the tissue comprising the airway wall and the parenchyma and *2*) the active force of ASM contraction.

We will first consider the passive forces (see Table 1 for definitions of symbols). When the ASM is relaxed, the tendency of the airway wall (which we consider to include the ASM itself) to recoil inward is balanced by the outward pull of the parenchyma. The parenchyma thus generates a transmural pressure (Ptm) across the airway wall. If the airway expands isotropically with pressure at the same rate as the parenchyma, then the parenchyma around the airway will remain undistorted as the lung is inflated. If we further neglect any axial pressure gradients along the airway, then Ptm will always equal Ptp. For this situation of isotropic parenchymal expansion to occur, the airway radius (*r*_{iso}) must increase with the one-third power of lung volume (20), which in the case of a constant lung elastance means that (1) where *r*_{isoTLC} is the value of *r*_{iso} at total lung capacity (TLC), and Ptp_{TLC} is Ptp at TLC.

In general, of course, the elastic properties of the airway wall are not perfectly matched to those of the parenchyma, in which case airway radius (*r*) and *r*_{iso} are not equal. This causes the parenchyma in the vicinity of the airway to be distorted, resulting in an additional component to the Ptm (ΔPtm), as indicated in Fig. 1. The relationship between ΔPtm and local parenchymal distortion caused by the airway depends on the shear modulus (μ) of the parenchyma. Lai-Fook (27) used a continuum mechanics approach to show that ΔPtm/2μ is a single-valued curvilinear function of *r*/*r*_{iso} [see Fig. 8 of Lai-Fook (27)]. Using a value for μ of 0.7Ptp (27), and multiplying by 1.4, allows ΔPtm/Ptp to be plotted against *r*/*r*_{iso} as a master curve from which ΔPtm can be determined for a given value of *r*/*r*_{iso} at any value of Ptp. This master curve is shown in Fig. 1*B* of Gunst et al. (20) and yields, by curve fitting to the original data of Lai-Fook (27), the approximate relation (2)

Combining *Eqs. 1* and *2* gives (3)

Therefore, (4)

The outward pull of the parenchyma on the airway wall is balanced by the inward recoil of the airway itself. This recoil arises from the elastic properties of the airway wall, which produce a passive hoop stress (*F*_{P}) when the wall is stretched past its point of elastic equilibrium. If we assume that *r* is substantially greater than the thickness of the airway wall (i.e., the airway behaves like a thin shell) and that the wall itself has negligible bending rigidity, then the relationship between *F*_{P} and the inward recoil pressure is given by the Laplace law for a cylinder. Thus: (5) Therefore, to calculate *F*_{P}, we need an expression that describes how *r* changes as a function of Ptm. We consider two possibilities. The first is that given by *Eq. 1* above, which describes an isotropically expanding lung, i.e., the airway expands at the same rate as the parenchyma when the lung is inflated. Equating *r* in *Eq. 4* to *r*_{iso} in *Eq. 1* gives the functional form for *F*_{P}, denoted *F*_{isoP}, that leads to isotropic expansion of the lung, namely (6)

When the ASM is stimulated to contract, it generates an additional active component, *F*_{A}, to the hoop stress in the airway wall, in which case *Eq. 5* becomes (7)

Solving *Eq. 7* for *F*_{A} and using *Eqs. 4* and *6* gives the expression for *F*_{A} corresponding to isotropic parenchymal expansion, denoted *F*_{isoA}, as (8)

Experimental evidence indicates, however, that the elastic properties of the parenchyma and airways are not perfectly matched. In general, therefore, *F*_{p} is not given by *Eq. 6* but instead assumes a functional form that causes the parenchyma around the airway to expand anisotropically as the lung is inflated. We will assume this functional form, denoted *F*_{anisoP}, to be that devised by Thorpe and Bates (45) to describe the canine airway pressure-radius data of Gunst et al. (18), namely (9) where *r*_{anisoTLC} is the airway radius achieved at TLC and *k* is a constant between 0 and 1. The airway wall described by *Eq. 9* can thus be viewed as having a circumference that is divided into two distinct regions. One region is completely inextensible and comprises a fraction *k* of the total circumference, whereas the remainder expands, as described by *Eq. 1*. Note that, even though *r*_{anisoTLC} is not, in general, equal to *r*_{isoTLC}, there will nevertheless be some value of Ptp for which the inflation of the airway precisely matches the inflation of the parenchyma (i.e., when ΔPtm is zero and Ptm equals Ptp). At this special value of Ptp, which we denote *P*_{0}, there will be no parenchymal distortion around the airway wall, so *r*_{iso} in *Eq. 1* and *r* in *Eq. 8* will be equal. Equating these two expression and solving for *r*_{anisoTLC} gives (10)

Combining *Eqs. 9* and *10* gives an expression for *F*_{P} under conditions of anisotropic parenchymal expansion, denoted *F*_{anisoP}; thus (11)

In analogy with *Eq. 8,* this gives a corresponding expression for the active ASM force, denoted *F*_{anisoA}, as (12)

Now, although the outward pull of the parenchyma is always balanced by the inward pressure generated by *F*_{A} and *F*_{P} via *Eq. 7*, this does not mean than *r* is constant. In fact, at any point in time, the ASM will be in the process of either lengthening or shortening depending on where *F*_{A} lies on the force-velocity relationship of the ASM. The length of the muscle is 2π*r*, so *r* will only be static if *F*_{A} is equal to the isometric force (*F*_{0}) of the muscle. In general, however, this will not be the case. The force-velocity relationship shows that *F*_{A} is greater than *F*_{0} when the muscle is lengthening (i.e., during eccentric contraction). Conversely, *F*_{A} is less than *F*_{0} when the muscle is shortening. We will assume that *F*_{A} is independent of ASM length and so depends only on the velocity of shortening according to the classic Hill equation (44), namely (13)

We used the above equations to calculate how *r* varies with time when the model is driven with a prescribed Ptp(*t*) signal, as follows. First, Ptp(*t*) was used in either *Eq. 8* or *Eq. 12* to determine *F*_{A} at each time step of 0.1 s. The result was substituted into *Eq. 13* to provide d*r*/d*t*, which was then used to determine *r* at the next time step using first-order Euler integration. This new value of *r* determines a new force balance condition yielding a new value for *F*_{A}, and so on. Finally, the calculated *r* profile was converted to a normalized airway resistance (Raw) by assuming Poiseuille flow in the airway and neglecting the effects of length change; thus (14)

#### Model fitting.

We tested the behavior of the model against previously reported time courses of Raw following intravenous methacholine injection in rats (7) and rabbits (38). The experimental Raw values were obtained by digitizing the published data at regularly spaced intervals from the printed pages, starting from when Raw began to increase above baseline. This yielded a set of data points for the rat and the rabbit to which the model was fit by least squares using a grid-search procedure (see appendix). The model was driven by a volume signal that varied sinusoidally above functional residual capacity so that Ptp in *Eq. 8* or *12* was given by (15) where *f* is ventilation frequency, *E*_{L} is lung elastance, and V_{FRC} is lung volume at FRC.

## RESULTS

Figure 2A shows experimentally measured time courses of Raw in rats following an intravenous bolus injection of methacholine when the lungs were held at constant volume. These volumes were established by setting the preconstriction positive end-expiratory pressure (PEEP) levels to 2, 4, and 6 cmH_{2}O, respectively (7). Also shown in Fig. 2*A* are the best-fit Raw time courses predicted by the model obtained at the same PEEP levels using *F*_{isoA} (*Eq. 8*), i.e., when the elastic properties of the airway wall were described by *Eq. 1*. Note the substantial systematic differences between the simulated curves and the corresponding experimental data points. The best-fit parameter values using *F*_{isoA} are listed in Table 2 along with the sum of the squared residuals (SSR) between the simulated and experimental data.

Figure 2*B* shows the Raw time courses obtained by Shen et al. (38) following an intravenous bolus of methacholine in rabbits. Here, the animals were ventilated, in each case at the same PEEP of 2 cmH_{2}O, but Vt was varied between 0, 5, and 10 ml/kg. (The error bars indicate the widths of the published curves caused by the flow oscillations used to estimate Raw. The smaller flow amplitudes used with the rats did not cause significant oscillations in the rat data used in Fig. 2*A*). Also shown in Fig. 2*B* are the best-fit model predictions of Raw obtained using Vt values of 0, 0.63, and 1.20 (expressed as fractions of FRC), which are similar in relative terms to those used to obtain the experimental data. The simulated curves show an inverse dependence of Raw on Vt as seen in the experimental data. However, as with the rat data in Fig. 2*A*, the simulated curves also exhibit marked deviations away from the experimental data. The best-fit parameter values and SSR are listed in Table 2.

We next repeated the above model-fitting exercise using *F*_{anisoA} (*Eq. 12*), i.e., when the elastic properties of the airway wall were described by *Eq. 9*. The results for the rat and rabbit data are shown in Figs. 3, *A* and *B*, respectively. In both cases, the simulated Raw curves fit the experimental data much better than when *F*_{isoA} was used (compare with Figs. 2, *A* and *B*); in the case of the rabbit data (Fig. 3*B*), the simulated curves intersect the error bars for most of the experimental data points. The best-fit values of the model parameters and the values of SSR obtained using *F*_{anisoA} are given in Table 2. The 5% significance value of the *F* ratio with 15 and 13 degrees of freedom, corresponding to the isotropic vs. the anisotropic model (17 independent data points minus 2 and 4 free model parameters, respectively) is ∼2.5. Because SSR obtained using *F*_{anisoA} was more than an order of magnitude smaller than that obtained using *F*_{isoA} in both species, a highly significantly improved model fit is provided by allowing the airway wall to have elastic properties that are different from those of the parenchyma.

Figure 4 shows the relationships between *r* and Ptm predicted by the two models of the airway wall, namely when *F*_{isoA} (*Eq. 8*) was used (identical for both rat and rabbit), and when *F*_{anisoA} (*Eq. 12*) was used. The *r* curves obtained using *F*_{anisoA} for both rat and rabbit are virtually identical and do not increase with Ptm nearly as quickly as does the *r* curve obtained using *F*_{isoA}. This shows that the rat and rabbit airways are considerably stiffer than the parenchyma.

## DISCUSSION

Although the main culprit responsible for the hyperresponsiveness of asthma is still controversial (11), an obvious possibility is that the ASM in asthmatic airways is somehow abnormal. Indeed, in vitro studies have shown that a key difference lies in the rate of shortening; asthmatic ASM shortens more rapidly than normal ASM (10, 32), and allergen-sensitized ASM contracts more rapidly than control ASM (26, 34). There is also evidence that velocity differences may play a role in the different responses of normal and asthmatic lungs to deep inspirations (25). Contraction velocity cannot be considered in isolation, however, because velocity is determined by the force-velocity relationship of ASM and the load against which the muscle must contract in situ (9). The central importance of load is demonstrated at the level of the whole organ following airway challenge by a marked inverse dependence of the rate of bronchoconstriction on lung volume (6, 7, 12, 41). It has even been postulated that the cyclic loads of airway-parenchymal interdependence during normal breathing may be sufficient to prevent normal ASM from proceeding to full contraction, whereas a pathologically rapid muscle may manage to narrow the airways excessively (16, 17, 42).

Theories such as these are intuitively compelling but difficult to evaluate unless their consequences are assessed quantitatively. This requires a mathematical model encapsulating the key events involved. Although the key events in ASM contraction are still a matter of some debate (3, 14, 15, 18, 19, 33, 35, 40, 42, 47), we recently showed that a simple model consisting of a force generator in series with a nonlinear spring can accurately simulate reported experimental force-length loops obtained under steady-state conditions in activated ASM strips (5). For the purposes of the present study, we used our tissue strip model in a circular configuration to represent the ASM that surrounds an airway. We also reasoned that the elastic load opposing ASM contraction provided by the outwardly tethering alveolar attachments to the airway wall (Fig. 1) would significantly outweight any load provided by the series elastic element within the ASM tissue itself. Therefore, in the present study, we omitted the series elastic element and merely incorporated a Hill-type force generator to model the ring of ASM. Anafi and Wilson (3) recently simulated the response of a constricted airway to a deep breath using a computational model similar to ours in terms of its mechanical linkage between airway and parenchyma. Their model uses a substantially more complicated, yet still entirely empirical, description of ASM contraction dynamics. We did not feel this extra complexity was warranted for the present application in view of the fact that we were recently able to accurately model the force-length behavior of an activated oscillating ASM strip using the relatively much simpler dynamics embodied in *Eq. 13* (5).

We first examined how the constricting airway model performed when simulating in vivo bronchoconstriction data for which the ASM was only in the process of shortening. The experimental data set we used for this purpose came from a previous study (7) in which we tracked the development of bronchoconstriction in rats following a bolus intravenous injection of methacholine while lung volume was held fixed. There were no cyclic changes in ASM length involved in this experiment apart from the small oscillations associated with the high-frequency measurement of Raw. We found that, provided the model takes airway wall stiffness explicitly into account, the model accurately reproduces the experimentally observed dependence of Ptp on the time course of bronchoconstriction (Fig. 3*A*). We next applied the model to a data set obtained in rabbits in which ASM length was sinusoidally oscillated with different amplitudes achieved by mechanically ventilating the animals with different Vt (38). Here, again, the model captures the essential features of the dependence of Raw on time and volume (Fig. 3*B*), despite the fact that the ASM switched continually between lengthening and shortening. Thus, even though the Hill equation was originally developed to describe force-velocity data obtained from quick-release experiments (23) in which muscle length changes in only one direction, we have found here that this equation can be used to give a good account of the dynamics of bronchoconstriction both when airway radius is shortening only (the rat data) and when it is cyclically shortening and lengthening (the rabbit data).

The above results thus suggest that the dynamics of bronchoconstriction can be understood to a large degree in terms of the relative stiffnesses of the parenchyma and airways, and the classical hyperbolic force-velocity relationship that was first described for skeletal muscle by Hill and that has more recently been applied to smooth muscle (9, 23, 44, 48). These conclusions, however, are based solely on model fits to data from whole animals. In all likelihood, the reality of the situation is not this simple, particularly as far as the smooth muscle is concerned. Numerous recent in vitro studies at the level of the single cell or tissue strip have shown that the dynamics of ASM contraction may be rather more complicated than the Hill equation would imply. There is strong evidence, for example, that the number of cross bridges taking part in ASM force generation at any point in time is highly dependent on the strain rate of the ASM (4, 22). This leads to the prediction of markedly hysteretic force-velocity behavior (15, 16, 28, 33, 35). Furthermore, the contractile apparatus in ASM appears to be able to adapt to the length at which the ASM is activated, due to plasticity in both the assembly of the contractile proteins and the points at which these proteins anchor to the cytoskeleton (19, 40). Nevertheless, we are still able to describe whole lung constriction dynamics using a simple tube of ASM obeying the Hill equation, perhaps because the effects of the many complexities mentioned above average out to produce an apparently much simpler ensemble behavior. It is also possible that the success of our model reflects the particular data sets we have considered here. Certainly, these data reflect typical experimental scenarios used to study bronchoconstriction dynamics in animals. On the other hand, this does not guarantee that the model will account for all conceivable types of experimental maneuvers designed to expose the transient nature of airway narrowing. Looking for situations in which our model fails will be important for future studies aiming, for example, to show that this simple model must be extended to include some of the complexities identified at the tissue level in vitro (18, 21, 40). In any case, it is clear that airway narrowing in vivo reflects a dynamic balance between the ability of the ASM to shorten itself and the ability of the opposing load to lengthen it, averaged over the period of interest. This may be, as Fredberg and coworkers suggest (15–17), a crucial determinant of bronchial hyperresponsiveness.

We must also consider the potential significance of other assumptions made in the model. Key among these is that force does not vary with ASM length. This not only simplifies the model by obviating the need to assign force-length properties to the ASM, but it also avoids the even more problematic issue of deciding how the length of ASM at FRC relates to the force-length curve. One might suspect, on teleological grounds, that ASM ought to be at optimum length for force generation at FRC. However, ASM is apparently able to change its optimal length depending on the length at which it is initially stimulated (19, 21, 24). This makes it difficult to know how the force-length curve of ASM should vary at different PEEP levels. Fortunately, the width of the peak in the ASM force-length curve is rather wide (44), so we felt that keeping maximal force independent of lung volume was a reasonable first approximation.

We also assumed that the activation of the ASM was immediate and uniform, so that when the model airway started to narrow it did so according to a force-velocity relationship that remained constant over the subsequent course of the narrowing. The situation is more complicated in reality because an intravenous injection of contractile agonist does not reach and activate all ASM in the lung at the same instant. Thus the experimental time courses of Raw shown in Figs. 2 and 3 likely reflect, to some degree, a distribution of agonist delivery rates throughout the lung. Indeed, this may explain the slight upward concavities seen in the very early stages of the experimental responses. Nevertheless, single airway sections in lung explants have been shown (52) to contract over a time scale similar to that of the bronchoconstriction data used in the present study, suggesting that the dynamics of distribution of an intravenous bolus of agonist probably have only a minor influence on the subsequent rate of bronchoconstriction. A related assumption in our model is that the smooth muscle activation caused by agonist injection does not decay with time. Such decay is inevitable and may be partly responsible for the plateaus in responses from the rabbits seen in Fig. 2*B* and 3*B*. In a previous study, for example, we modeled this decay in dogs using interconnected compartments to account for the dynamics of uptake and clearance of agonist (29).

Another major simplification in our model is that we represent the airway tree of the lung in terms of a single equivalent airway. The experimental measurements of Raw we tested our model against in the present study (7, 38) were all obtained using single-frequency oscillations of flow at frequencies significantly above that of normal respiration, and as such represent the resistance of the airway tree with little contribution from the tissues (6, 7). The amplitudes of the oscillations were also small, justifying our use of the Poiseuille relationship (*Eq. 14*) to calculate Raw. Nevertheless, modeling in vivo bronchoconstriction in an entire lung in terms of what happens to a single airway neglects the fact that airway narrowing is invariably heterogeneous (30, 31, 45). Thus, although our model describes the data rather well using putative mechanisms that reflect ASM dynamics at the level of the isolated muscle strip, we have to bear in mind the possibility that whole lung behavior may also reflect a distribution of small-scale events such as airway closure or self-organized clustering of ventilation defects, as recently modeled by Venegas et al. (49).

We also assumed that the parenchyma surrounding this single airway is purely elastic and obeys the pressure-distortion relationship derived from the work of Lai-Fook (27) expressed in *Eq. 2*. This expression accounts for the nonlinear relationship between parenchymal distortion and recoil pressure but neglects viscoelasticity and assumes circular symmetry about the center of the airway. In reality, the parenchyma is viscoelastic, so the rate at which it is distorted during airway contraction will have an influence on its apparent stiffness. Also, parenchymal distortion around a contracting airway can be highly asymmetric as a result of heterogeneities in the structure of the nearby parenchyma (2). The single airway in our model thus does not correspond to any particular anatomical structure in the lung but rather acts as a vehicle for capturing the essential elements involved in airway narrowing.

An interesting consequence of fitting our model to experimental data is that we obtained an estimate of the pressure-radius relationship for the airway wall in vivo. Of course, this relies on the assumption that the elastic properties of the airway wall are described to a sufficient degree of accuracy by *Eq. 9*. This expression was developed by Thorpe and Bates (45) and is based on the canine data of Gunst et al. (18) and is convenient in that it is defined by only two free parameters, *P*_{0} and *k*. Even so, this expression is also likely useful for application to other species. In particular, Shen et al. (39) measured diameter-pressure relationships for airways of various generations in excised rabbit lungs and found qualitatively similar curves to those we estimated for rat and rabbit shown in Fig. 4. They also found that as Ptp varied between 0 and 20 cmH_{2}O, airway diameter changed by between ∼15 and 45% of its value at 20 cmH_{2}O depending on airway generation number. The corresponding change for the rat and rabbit curves estimated in the present study (Fig. 4) is ∼30%, suggesting that our curves represent a reasonable average for the airways of the rabbit. Sera et al. (37), on the other hand, measured airway distensibility in rats using microcomputed tomography and found, not surprisingly, that specific compliance varied widely throughout the lung and tended to be greater as airway size decreased. In any case, the functional form of *Eq. 9* merely implies that the circumference of the single equivalent airway in our model is composed of two components, one that is perfectly rigid and one that expands with the cube root of volume similar to *Eq. 1*, with the parameter *k* determining the relative roles of the two components. This is a simple but very general model of wall stiffness that does not have special relevance to any particular species. Furthermore, inclusion of this expression in the model gives a very significantly improved model fit compared with when the airway wall is elastically matched to the parenchyma (Fig. 3*A* vs. 2*A* and Fig. 3*B* vs. 2*B*), suggesting that *Eq. 9* captures an important aspect of airway wall behavior. Another assumption was that the airway can be treated like a thin shell, allowing the Laplace law to be used to relate *r* to Ptm through *Eq. 5*. Although this may be reasonable in a normal lung, it may be less so in allergic asthma when inflammatory thickening of the airway wall can be substantial (51). Pathologically remodeled airways might also have increased bending rigidity, also contrary to the assumptions implicit in *Eq. 5*. Nevertheless, the model based on *F*_{anisoA} may still prove useful for the study of the mechanical consequences of airway remodeling in conditions where wall stiffness may be expected to change.

Another consequence of fitting our model to experimental data is that we estimate values for the parameters *a*, *b*, and *F*_{0} that characterize the contractility of ASM via *Eq. 13*. These parameters cannot be evaluated in vivo directly, but they have been determined for isolated strips of ASM in vitro (23, 36, 43). We therefore need to ask whether the values we obtained for *a*, *b*, and *F*_{0} (Table 2) are compatible with those found in vitro. We forced this compatibility to a certain degree (see appendix) by requiring that *a* = 0.25*F*_{0}, as has been reported for rat ASM (9). Relating the values of *a*, *b*, and *F*_{0} to the known biophysics of ASM any further becomes problematic, however, because the values of *b* (which has units of velocity) and *a* and *F*_{0} (which have units of force) are functions of scale. For example, doubling the cross section of the muscle will double *F*_{0}, whereas doubling its length will double *b*. Another problem is that the values of *a*, *b*, and *F*_{0} reported in the literature come from maximally stimulated strips of ASM, whereas the rat and rabbit data we used in the present study were obtained using by submaximal doses of methacholine. Indeed, the rabbits (38) received a significantly smaller dose than the rats (7), which explains why the values of *F*_{0} from the anisotropic model estimated for the rabbits are smaller than those for the rats (Table 2). Nevertheless, we can proceed a little further by noting that in vitro studies estimate the unloaded contraction velocity, *bF*_{0}/*a*, of maximally stimulated ASM (43) to be ∼0.25 *l*_{0}/s, where *l*_{0} is initial muscle length. *l*_{0} in absolute units is thus equal to 4*bF*_{0}/*a* calculated using nonnormalized values for *a*, *b*, and *F*_{0}. When stimulation is submaximal, *bF*_{0}/*a* is <0.25 *l*_{0}/s, and *l*_{0} is correspondingly greater than 4*bF*_{0}/*a*. This suggests that 4*bF*_{0}/*a* calculated using the parameter values in Table 2 provides a lower bound on *l*_{0}, which for our model is equal to the circumference of the airway in Fig. 1. Accordingly, we used the values of *a*, *b*, and *F*_{0} for the anisotropic model in Table 2 to obtain lower bounds for *r*_{p} of 1.2 mm for the rat and 5.8 mm for the rabbit. Because we do not know how far from maximal stimulation the rats and rabbits were, we have no idea how much these values for *r*_{p} underestimate the radii of single airways that serve to represent the entire airway tree in each species. We can say, however, that at least they seem to be of the correct order of magnitude.

In conclusion, we have developed a model to account for the time course of bronchconstriction in terms of a single airway narrowing against the opposing elastic loads presented by the stiffness of the airway wall and the outward tethering forces of the parenchyma. Despite the simplistic structure of this model, it is able to simulate the major features seen in Raw following a bolus injection of methacholine as a function of both PEEP and Vt. Our results indicate that much of the dynamic nature of bronchoconstriction in vivo can be understood in terms of the classic Hill force-velocity relationship operating against the elastic loads presented by the surrounding lung parenchymal attachments and an airway wall that is stiffer than the parenchyma.

## APPENDIX

The experimental data sets used to test the model in the present study each consisted of three separate time courses for Raw. These time courses were obtained in the rat at three different levels of PEEP, whereas in the rabbit they were obtained by ventilating the animal with three different Vt. We fit two different models to each of these two data sets. The first is the isotropic model that incorporates *F*_{isoA} (*Eq. 8*) to describe the elastic properties of the airway wall and is characterized by the three parameters of *Eq. 13*, namely *F*_{0}, *a*, and *b*. The second is the anisotropic model that incorporates *F*_{anisoA} (*Eq. 12*) to describe the elastic properties of the airway wall and introduces the two additional parameters *k* and *P*_{0} through *Eq. 10*. We fit these two models to each set of Raw time courses by finding values for the parameter sets {*F*_{0}, *a*, *b*} and {*F*_{0}, *a*, *b*, *k P*_{0}}, respectively, that minimized the total SSR between the three experimental time courses and the predicted values of Raw at corresponding time points. The model fitting was achieved using a grid-search procedure in which the fits provided by the isotropic and anisotropic models were systematically examined over a grid with axes corresponding to values for each of the free parameters in each model. This did not include the parameter *a* in *Eq. 13* because we applied the constraint *a* = 0.25*F*_{0} to match in vitro measurements on ASM from rats (9).

Because *a* is determined by *F*_{0}, the isotropic model has only two free parameters, *F*_{0} and *b*. To fit this model to the data, initial search ranges for these two parameters were determined by experimenting with various ranges within which the eventual best-fit values lay. Suitable ranges were found to be 1–50 (arbitrary force units) for *F*_{0} and 0.001–4.0 (arbitrary velocity units) for *b*. Test values for each parameter were chosen at the borders of *N* equally spaced intervals along its range. This defined *N*^{2} sets of parameter values located at points on a grid in two-dimensional parameter space. We experimented with various values of *N* and established that *N* = 6 lead to consistent best-fit values for both *F*_{0} and *b*. To simulate the rat data with each set of parameter values, the model was driven via *Eq. 15* using *V*_{FRC} values of 400, 800, and 1,200 ml and Vt = 0. To simulate the rabbit data, *V*_{FRC} = 400 ml with Vt = 0, 250, and 500 ml. The remaining model parameters were assigned the values Ptp_{TLC} = 30 cmH_{2}O, *r*_{TLC} = 1 cm, and *E*_{L} = 5 cmH_{2}O. These values correspond to a normal lung in relative terms but represent an arbitrary length scale. To adjust for this arbitrary choice of scale, each simulated set of Raw curves was divided by the initial experimental value of Raw from either the PEEP = 2 cmH_{2}O data (for the rat) or the Vt = 0 data (for the rabbit). The value of SSR between each set of scaled simulated Raw curves and the corresponding experimental data points was then calculated, and the set of model parameter values yielding the minimum value of SSR was taken as the best-fit set. This best-fit set was then improved on by repeating the grid-search procedure 12 times, each time with the new grid being centered on the best-fit values from the previous search but the dimensions of the grid being 1/2 that of the previous grid, with the proviso that no parameter was allowed to become negative. The final spacing between consecutive values on the parameter grid was <1% of the best-fit values for each parameter.

The anisotropic model has an additional two parameters, *k* and *P*_{0}, accounting for airway wall stiffness. The parameter *k* must lie between 0 and 1, so we defined its initial search range to be 0.01–0.99. We considered that the parameter *P*_{0} should lie somewhere below Ptp at TLC, so we defined its initial search range to be 1–30 cmH_{2}O. We then performed a two-dimensional grid search with *N* = 6 for these two parameters. At each of these points, a grid search on *F*_{0} and *b* was performed as described above for the isotropic model. The refinement of the *k*-*P*_{0} grid was performed 12 times, as with the first model, by halving the grid dimensions at each pass. Fitting the anisotropic parameter model by grid search was not, however, entirely straight-forward. Preliminary investigations with various combinations of values for *F*_{0}, *b*, *k*, and *P*_{0} in both species showed that the five-dimensional hyper-surface defined by SSR had a complicated shape without a single clearly defined minimum (Fig. 5). Consequently, the best-fit values obtained for *F*_{0}, *b*, *k*, and *P*_{0} were somewhat sensitive both to the values of *N* used to define the grid points and on the number of times the grid dimensions were refined. Experimentation with various combinations of these quantities indicated that the values we used gave parameter estimates close to the global minimum, which is supported by visual assessment of goodness-of-fit in Fig. 3.

Something of particular note that we found in fitting the anisotropic model to the data was that *F*_{0} and *P*_{0} were able to move substantially in opposite directions without affecting the value of SSR very much. This can be appreciated in Fig. 5, which shows the surfaces defined by SSR for various values of *k* and *P*_{0}. In both surfaces, there is a trough positioned around *k* = 0.6 along which *P*_{0} can take any value over its entire range without having much effect on SSR. Therefore, we cannot place too much stock in the absolute values of either *P*_{0} or *F*_{0}. Nevertheless, *k* was relatively insensitive to the details of the grid search, as shown by the narrowness of the troughs in Fig. 5. Fortunately, it is the value of *k* that is most relevant for the stiffness of the airway in terms of how much *r* increases for a given increase in Ptp, as shown in Fig. 4. Accordingly, we feel that our model-fitting procedure is able to robustly estimate airway wall stiffness from experimental data such as that used in the present paper.

## GRANTS

This work was supported by National Institutes of Health Grants R01 HL-67273, R01 HL-75593, and P20 RR-15557, and by the Fonds de le Recherches en Sante du Quebec.

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