## Abstract

During methacholine challenge tests of airway responsiveness, it is invariably assumed that the administered dose of agonist is accurately reflected in the dose that eventually reaches the airway smooth muscle (ASM). However, agonist must traverse a variety of tissue obstacles to reach the ASM, during which the agonist is subjected to both enzymatic breakdown and removal by the bronchial and pulmonary circulations. This raises the possibility that a significant fraction of the deposited agonist may never actually make it to the ASM. To understand the nature of this effect, we measured the time course of changes in airway resistance elicited by various durations of methacholine aerosol in mice. We fit to these data a computational model of a dynamically contracting airway responding to agonist that diffuses through an airway compartment, thereby obtaining rate constants that reflect the diffusive barrier to methacholine. We found that these barriers can contribute significantly to the time course of airway narrowing, raising the important possibility that alterations in the diffusive barrier presented by the airway wall may play a role in pathologically altered airway responsiveness.

- airway resistance
- hyperresponsiveness
- computational model
- mouse model

a variety of disparate mechanisms have been invoked to explain measurements of airway hyperresponsiveness (AHR) (2, 7, 11). Virtually all, however, make the assumption that the administered dose of agonist is accurately reflected in the dose that eventually reaches the airway smooth muscle (ASM). In fact, this is not the case. It is well appreciated, for example, that the fraction of inhaled agonist that is deposited on the airway walls and the sites of deposition are markedly influenced by the pattern of respiratory flow, droplet size, and the geometry of the airway tree (20, 22). These can both change substantially as bronchoconstriction proceeds. Indeed, it has been suggested that the response elicited by aerosol challenge is self-limited due to the fact that, when airways become severely narrowed, they prevent further accumulation of agonist (10).

Reaching the airway wall, however, is not the only challenging leg of an agonist's journey to the ASM. Once impaction has occurred, the agonist must traverse layers of fluid and mucus, epithelium, basement membrane, connective tissue, and interstitium to reach the ASM and elicit a response. During this diffusive transport phase, the agonist is subjected to both enzymatic breakdown and removal by the bronchial and pulmonary circulations (15, 16, 26, 27). This raises the possibility that a significant fraction of the deposited agonist may never actually make it to the ASM, which, in turn, suggests that changes in the diffusive barrier of the airway wall will affect airway responsiveness. Indeed, our laboratory has previously reported that the time to peak response in aerosol-challenged mice is increased when the mice are allergically inflamed (24, 25) and decreased when their epithelial barrier function is compromised by treatment with cationic protein (4, 8).

The diffusive transport of agonist across the airway wall thus has substantial implications for linking the phenomenon of AHR to its underlying mechanisms in animal models. Accordingly, the goal of the present study was to gain a quantitative understanding of how the kinetics of methacholine transport affect airway responsiveness in mice. We pursued this goal by measuring the time course of bronchoconstriction elicited by various durations of methacholine aerosol and interpreted the results in terms of a computational model of a dynamically contracting airway responding to agonist that diffuses through an airway compartment. Fitting this model to experimental data provides values of intercompartmental rate constants that reflect diffusive barriers to methacholine, as well as an estimate of the fraction of the originally delivered dose of methacholine that eventually reaches the ASM.

## METHODS

These studies were approved by the Institutional Animal Care and Use Committee of the University of Vermont and were in accord with the animal welfare act.

#### Animal preparation.

We studied two groups of BALB/c mice at 8–12 wk of age and 18–22 g (Jackson Laboratories, Bar Harbor, ME) that were sensitized with two intraperitoneal injections of ovalbumin (100 μg) in alum (2.25 mg), given on *days 1* and *14*, after first having been acclimatized in our animal facility for 3 days. One group (the control group) was then exposed to an aerosol of saline for 30 min on *days 21*, *22*, and *23*. The other group (the OVA group) was challenged by exposure to an aerosol of 1% ovalbumin for 30 min on *days 21*, *22*, and *23*. The experimental protocol (see below) was performed on both groups on *day 25*.

#### Experimental protocol.

The animals were anesthetized (pentobarbital sodium, 90 mg/kg ip) and tracheostomized, and an 18-gauge cannula was tied into the trachea. The animals were then attached to a Flexivent mechanical ventilator (Scireq, Montreal, Quebec, Canada) and ventilated with a piston displacement of 0.25 ml (volume delivered to the lungs ∼0.20 ml after gas compression in the ventilator cylinder is accounted for) at a frequency of 200 breaths/min against a positive end-expiratory pressure (PEEP) of 3 cmH_{2}O. An intraperitoneal injection of pancuronium bromide (0.8 μg/kg) was given for paralysis after which the ECG was monitored as an autonomic indicator of appropriate depth of anesthesia.

Each mouse was given a deep lung inflation to a pressure limit of 25 cmH_{2}O to recruit the lungs and standardize volume history. After approximately 1 min more of mechanical ventilation, the animals were allowed 1 s for expiration and then respiratory input impedance (Zrs) was determined with the application of a 2-s broadband volume perturbation signal. The volume perturbation contained 13 mutual prime frequencies from 1 to 19.625 Hz with amplitudes that scaled approximately inversely with frequency and had a peak-peak piston displacement of 0.17 ml. The animals were then exposed in groups to either two, four, six, or eight consecutive methacholine challenge epochs (*n* = 6, 8, 6, and 6, respectively, for the control group, and *n* = 6, 7, 6, and 8, respectively, for the OVA group). Each challenge epoch began with 10-s exposure to a methacholine aerosol (3.15 mg/ml), followed immediately by the measurement of Zrs. Note that the events involved in the measurements of Zrs included turning off the aerosolization of methacholine after it had been running continuously for 10 s, allowing exhalation to functional residual capacity for a period of 1 s, delivering a 2-s volume perturbation, and then finally turning the aerosolization back on again. These events took ∼7 s, so a single complete epoch of aerosolization and Zrs measurement took ∼17 s. Thus, for example, an animal receiving two methacholine challenge epochs would be exposed to aerosol over two consecutive 10-s periods separated by the ∼7 s involved in making the intervening measurement of Zrs.

Once the challenge epochs had been completed, the regular measurements of Zrs continued for the next 3 min. In this way, we determined how Zrs changed both while the sequence of methacholine aerosolizations were being delivered and then subsequently. Note that it was necessary to run the animals in separate cohorts rather than having each animal be subjected to each duration of methacholine aerosolization, because we found that, although the animals responded reasonably uniformly to a first challenge, their responses to subsequent challenges were much more variable, possibly as a result of secretions produced by the first challenge.

The constant-phase model of impedance (14) was fit to each measurement of Zrs to yield values for the Newtonian resistance of the respiratory system, which we can take to be a measure of overall airway tree resistance (Raw) (23). The constant-phase model also provides measures of tissue viscance (G) and stiffness (H), which we know from previous work from our laboratory (24) reflect events taking place in the lung periphery. In particular, increases in H signify closure of small airways due to liquid bridging across the lumen (19). Although these important physiological events do not relate directly to the general narrowing of the conducting airways, they do provide important information about the extent to which the lung behaves homogeneously (3).

#### Computational modeling.

To quantitatively interpret the data collected in the above experiments, we developed a computational model that we believe embodies the key physical processes involved in the observed dynamics of Raw following methacholine aerosol exposure. The active component of the model consists of a circular elastic airway embedded in uniform elastic parenchyma, as our laboratory has previously described (6, 12). The airway wall is considered to contain a ring of ASM that contracts, when stimulated, according to the classic Hill hyperbolic force-velocity relationship (13). As the ASM in the wall shortens, opposing elastic stresses are generated in the wall itself and in the attached parenchyma that becomes increasingly distorted as the airway narrows. These opposing forces determine the load on the ASM, which, in turn, determines its contraction velocity through the Hill equation. The result is a time course of airway narrowing that starts from the relaxed configuration and proceeds until the forces of ASM contraction are no longer able to overcome the steadily increasing opposing elastic forces. This model contains the following parameters: ASM isometric force (F_{0}), the parameters *a* and *b* in the Hill equation that define the hyperbolic shape of the ASM force-velocity relationship (13), a parameter *k* that determines the stiffness of the airway wall relative to the intrinsic stiffness of the surrounding parenchyma (0 < *k* < 1), and a parameter P_{0} that specifies the transmural pressure across the airway wall at which the airway and parenchymal inflation states are matched (i.e., where parenchyma attached to the wall is not distorted). We have found that this model can be accurately fit to experimental measurements of Raw vs. time following intravenous injection of methacholine at three different levels of PEEP (1, 3, and 6 cmH_{2}O) by adjusting only three free parameters: F_{0}, *b*, and *k* (6, 12). The remaining two parameters can be prespecified; the model fit is insensitive to the value of P_{0} (6), so we set this parameter equal to 3 cmH_{2}O so as to place it in the middle of the PEEP range that we investigated, and it has been shown that *a* is numerically equal to one-fourth of the value of F_{0} (9).

This airway model as it stands, however, is not sufficient to account for the aerosol exposure data of the present study, because it assumes that the ASM becomes maximally activated instantaneously as soon as methacholine is administered. As such, it was suitable for simulating the bronchoconstriction that results from a bolus intravenous injection of agonist (5, 6, 12), but it cannot account for the extended duration of agonist delivery that occurs during aerosolization. For this, we need to incorporate a description of the progressive accumulation of agonist at the epithelial surface and its subsequent diffusion through the tissues that separate the airway lumen from the underlying ASM. Our laboratory has shown previously in dogs (18) that the slower kinetics of bronchoconstriction following aerosol delivery relative to intravenous injection can be accounted for by assuming that the airway wall acts like a volume of distribution from which agonist diffuses at a rate proportional to its concentration within the compartment. Therefore, in the present study, we combine this compartment modeling approach with our dynamically contracting airway model to account for the kinetics of bronchoconstriction, resulting from aerosol delivery of methacholine in mice.

The model of agonist transport consists of airway wall and tissue compartments, as shown in Fig. 1. The deposition of agonist into the airway wall compartment by aerosolization is represented by the function *a*(*t*), which has a constant value for the duration of aerosolization. This causes the concentration, C_{w}(*t*), of agonist in the airway wall compartment to rise as agonist accumulates within it. At the same time, agonist is exchanged between the airway wall and tissue compartments by passive diffusion, driven by the bidirectional rate constant, *k*_{dif}. When aerosolization is terminated, *a*(*t*) becomes immediately zero, while diffusive exchange between the two compartments continues. The concentration of agonist, C_{t}(*t*), in the tissue compartment is determined by the balance between its rate of diffusion in from the airway wall compartment and its rate of diffusion out both to the airway wall compartment and to a sink. The latter accounts for the permanent removal of agonist, which is presumed to occur through a combination of enzymatic digestion and clearance by the local circulation. It is not important for the model to distinguish between these two removal processes, although, if the situation in mice is similar to that in dogs (18), we suspect that blood flow is the major sink mechanism. In any case, the rate constant of the sink is *k*_{sink}. C_{t}(*t*) is taken to be concentration of agonist experienced by the ASM and is, therefore, proportional to the level of ASM activation, which itself is determined by an activation factor (γ) that scales the values of F_{0}, *a*, and *b* in the ASM force-velocity relationship. The complete set of equations describing the entire model is given in the appendix.

#### Model fitting.

We fit the computational model to the complete data sets obtained from each of the two groups of mice studied (the control group and the OVA group). That is, the values of the model parameters were determined that gave the best (least squares) fit to the data of Raw vs. time for all four durations of aerosolization simultaneously. To keep the free parameters in the model to a manageable number, we set the values of the Hill equation equal to those found in a previous study for normal BALB/c mice (12) (F_{0} = 2.7 cm/cmH_{2}O, *a* = F_{0}/4, *b* = 0.03 cm/s). Similarly, we set the airway wall stiffness parameter (*k*) equal to 0.8. This left three free parameters, namely *k*_{dif}, *k*_{sink}, and γ. We determined the best fit values of these parameters using a grid search procedure in which initial wide ranges were specified for each parameter within which seven equally spaced values were chosen. The model was fit to the data using every combination of the three parameter values, and the combination giving the smallest value of the root mean squared residual (RMSR) was retained. The search range was then reduced to one-half its previous size, centered about the best fit values, and the search was repeated. This process was repeated five times. The model was forced to have a baseline Raw value of 0.3 cmH_{2}O·s·ml^{−1} to match the average experimental value and was fit only to the early portions of each data set over which Raw continued to increase, because, once the peak had been reached, the data descended to an intermediate plateau level that likely reflects phenomena not included in the computational model, such as peripheral airway closure (19, 24).

## RESULTS

Figure 2 shows the individual data sets obtained from each of the mice studied, to demonstrate that, although there was some variation in the magnitude of the responses obtained with each duration of methacholine aerosolization, the timing of the peak in the response was similar in all cases. (Some of the data points are missing in the figure for the two longest aerosol durations because the large changes in Zrs caused its fit by the constant-phase model to be poor, but this applied to <10% of the data.) The baseline values of Raw averaged across all mice were 0.347 (SD 0.092) cmH_{2}O·s·ml^{−1} for the control group and 0.383 (SD 0.092) cmH_{2}O·s·ml^{−1} for the OVA group. These baseline values were not significantly different between the two groups (unpaired *t*-test, *P* > 0.05). Raw following methacholine aerosolization was significantly different between the two groups (ANOVA, *P* < 0.05) for aerosolization durations of 34, 68, and 136 s, but not for 102 s.

Figure 3 shows Raw, G, and H from all groups studied. Of note, all three parameters followed qualitatively similar time courses, with the relative responses in G and H being as great as those in Raw.

Figure 4 shows the time courses of Raw measured in the control group mice during and following 34, 68, 102, and 136 s of methacholine aerosolization. Note that the total dose of methacholine delivered to the animals increases with the duration of aerosolization, so the peak responses in Raw increase accordingly with duration (Fig. 2). Also shown are the fits provided by the computational model, obtained by adjusting the values of *k*_{dif}, *k*_{sink}, and γ through the grid search procedure described above in methods. Although the model fits show systematic deviations on either side of the data, they nevertheless recapitulate, in an overall sense, the way that Raw increases over time with a response amplitude that increases with the duration of aerosolization. Importantly, the fits also mimic the way that Raw continues to increase for a while after termination of methacholine delivery.

We also found that the quality of the model fit was extremely insensitive to the values of *k*_{sink} and γ, but it was very sensitive to the ratio of these two parameters. Specifically, these two parameters could each be varied about their best fit values by an order of magnitude without decreasing the quality of the model fit (i.e., increasing the value of RMSR) by more than ∼1.5%, provided the ratio *k*_{sink}/γ varied by <10%. Similarly, *k*_{dif} could only vary by ∼18% around its best fit value without increasing RMSR by 1.5%. In other words, the model effectively has only two independent parameters, *k*_{dif} and the ratio *k*_{sink}/γ, and so we report the values of these two quantities in Table 1, along with the RMSR value for the fit.

Figure 5 shows corresponding plots for the OVA group, while the values of the best fit model parameters and the RMSR are again listed in Table 1. Notably, the value of *k*_{dif} in the OVA group is about one-third its value in the control group, suggesting that the diffusivity of methacholine across the airway wall is substantially decreased in the OVA group. The ratio *k*_{sink}/γ is also about three times smaller in the OVA group compared with the control group. The computational model also provides, through the integration of *Eqs. A4* and *A5* using the parameter values in Table 1, a means for estimating the relative contributions of methacholine that reach the ASM in the control vs. the allergically inflamed animals. Due to the smaller value of *k*_{sink}/γ in the latter group, these equations predict that methacholine concentrations reach as much as 64% higher in the inflamed animals over the course of the longest duration of aerosol delivery.

The above results from the OVA group would seem to suggest that allergic inflammation had the effect in these mice of increasing the diffusion barrier presented by the airway wall and of making the ASM more sensitive to methacholine. However, our laboratory has reported previous evidence (19, 24) that the ASM in allergically inflamed BALB/c mice behaves normally, and that the hyperresponsiveness of these animals is due almost entirely to a physically thickened epithelium that geometrically amplifies a normal degree of ASM shortening. To test whether the same might apply in the present study, we refit the computational model to the data from the OVA group, but this time the model was given a lining inside its airway to simulate the effects of a thickened epithelium (see appendix). We found, by trial and error, when this lining had a thickness of 2% of the relaxed airway diameter, we obtained best fit values for *k*_{dif} and *k*_{sink}/γ that were very similar to those found for the control group (Table 1), with minimal degradation in the quality of the model fit (Fig. 6). Furthermore, with these parameter values, the prediction of the time course of C_{t}(*t*) made by the computational model (*Eqs. A4* and *A5*) indicates that the methacholine concentration reaching the smooth muscle is now only as much as 5% greater than in the control animals.

## DISCUSSION

Our study was motivated by the notion that the finite time required for an agonist to diffuse across the airway wall plays a significant role in the kinetics of bronchoconstriction. Evidence for this notion in mice comes from the observation that allergic inflammation seems to retard the peak in Raw produced by aerosol methacholine challenge (24, 25), whereas intratracheal instillation of cationic protein causes the peak to be achieved earlier than in control animals (8). We have also shown that bronchoconstriction kinetics in dogs are more drawn out in time when agonist is administered by aerosol vs. by intravenous injection (18). In the latter case, we applied a compartmental analysis to the time courses of bronchoconstriction obtained by aerosol vs. injection, which led us to conclude that the time constant of transport across the airway wall is in the order of 1 min. In the present study, our goal was to apply a related model analysis to the time course of bronchoconstriction in mice, and in particular to determine whether this analysis indicates that the airway remodeling associated with allergic inflammation causes an increase in the diffusive barrier of the airway wall, as our laboratory has previously suggested (24, 25). The importance of such a diffusive barrier stems, in large part from the fact that, while methacholine is in the process of diffusing across the airway wall, it is also being subject to both enzymatic degradation and removal by the local circulation (15, 16). This means that not all of the agonist that is deposited on the airway mucosa makes it through to the ASM. Consequently, any process that affects diffusion time will also alter the fraction of agonist that eventually stimulates the ASM to contract.

We used a simple yet novel experimental protocol in the present study, namely to measure the kinetics of Raw both during and after administration of methacholine aerosol into the airways. We hypothesized that continued passage of methacholine through to the underlying tissues from its site of initial buildup on the surface of the airway wall would cause it to have a concentration-time profile at the ASM that is delayed relative to its delivery into the airways. In other words, the peak in Raw ought to be delayed relative to the delivery of the aerosol, and the magnitude of this delay should be indicative of the characteristic diffusion time across the airway wall. This is clearly borne out by our data; Raw exhibits a continuing upward trajectory for ∼30 s following cessation of methacholine aerosol delivery (Figs. 2⇑⇑⇑–6). Furthermore, the inverse of the estimated rate constant of transfer between the airway wall and tissue compartments (i.e., *k*_{dif}, Table 1) equals the time constant of diffusion of methacholine across the airway wall and has a value of ∼80 s. This is close to the rate constant we determined previously for dogs (18), which would seem to suggest that the diffusive properties of the airways in these two species are similar, a curious finding given their great differences in lung size. We can only speculate on why this similarity should exist, but it is perhaps expected because rates of clearance of inert tracer from the lungs have been shown to be of the same order over a wide range of species size (28).

The fact that we see functional evidence of methacholine transport kinetics across the airway wall is itself hardly a surprise. What is intriguing, however, is whether these kinetics are altered in diseases that involve airway remodeling. To address this question, we repeated our analysis in a standard mouse model of allergic inflammation, which exhibits markers of remodeling, at least of a transient nature, such as a thickened epithelium and increased airway mucus (24). Our initial investigation indicated that this question is answered in the affirmative; when we fit the model to the data from the OVA group, we found that *k*_{dif} was reduced to ∼40% of control (Table 1). We also found that the ratio *k*_{sink}/γ was roughly one-third that of control, which would suggest that the activation of the ASM in the OVA group was elevated. This presents a problem, however, because our laboratory found clear evidence in a previous study (24) that the activation of the ASM in allergically inflamed BALB/c mice treated with methacholine aerosol is of a normal level, and that the AHR in this preparation is due to a physically thickened airway wall, which causes a geometric amplification of the degree of luminal narrowing and small airway closure (19, 24).

We, therefore, added a lining to the airway wall in the computational model so that the same degree of increase in Raw could be achieved with a smaller level of ASM activation. Specifically, a lining thickness of 2% of the relaxed airway diameter was sufficient to increase the ratio *k*_{sink}/γ back up to the level found in the control group. At the same time, however, this resulted in a value of *k*_{dif} that was close to that of the control group (Table 1). It is difficult to know how to interpret a 2% lining thickness in anatomic terms in a computational model that includes only a single airway, because, in reality, a thickening of the wall could apply to every branch of the airway tree. When we modeled this effect in a previous study using an anatomically based model of the airway tree, we found that a lining of 18 μm in every airway branch was sufficient to account for the wall thickening effect on airway responsiveness (24). This would correspond to a 2% increase in wall thickness in an airway of 1-mm diameter. Interestingly, this is about the size of the largest airways in a BALB/c mouse (21), which perhaps can be taken to indicate some level of consistency. This being the case, our results concur with those of our laboratory's previous study (24), which led to the conclusion that the hyperresponsiveness of the acutely allergically inflamed BALB/c mouse can be explained virtually entirely as being due to a thickening of the airway epithelium. On the basis of the present study, we can add the further conclusion that these changes in the airway wall are not accompanied by detectable changes in the diffusive barrier it presents to methacholine.

Thus, although our experimental data make it clear that the diffusive delay of methacholine delivery to the ASM is a significant component of the kinetics of bronchoconstriction, it is also clear that one cannot determine whether there are any changes in the diffusive delay without taking into account any changes in airway geometry. Nevertheless, it is also clear that changes in the diffusive transport of agonist have the potential to contribute to aberrations in airway responsiveness and might play a significant role in diseases involving significant airway remodeling. Diffusive transport thus joins the other recognized classes of biophysical mechanism that potentially play a role in AHR (2, 11).

Our model-fitting approach to analyzing the time course of bronchoconstriction has led us to a number of physiological conclusions. It is important to bear in mind, however, that these conclusions only hold to the extent that our computational model captures the essential features of reality, and on this score there are a number of important caveats. In particular, there is the structural simplicity of the model to consider. We represent the lung as a single contracting airway, while the diffusive pathway of agonist to the ASM is represented as only two interacting compartments (Fig. 1). These are the simplest representations we could devise that are able to capture the essential functional attributes of a dynamically contracting airway separated from its agonist by a semipermeable barrier, but they are still severely simplifying approximations of the relevant structures in a real lung. We also assumed that the lung behaves linearly in its response to methacholine in the sense that rate of delivery of methacholine to the surface of the airway wall was constant during the period of aerosolization, which does not account for the possibility that delivery might become reduced as the airway narrows (1). In fact, the concurrent increases in G and H (Fig. 3) show that significant numbers of small airways must have become closed during the methacholine response, particularly in the allergically inflamed animals, as our laboratory has shown previously (19, 24). Nevertheless, we neglect these effects in our computational model. Such simplifications are necessary in an inverse model if its parameters are to be uniquely identified from a limited data set. This necessary model simplicity likely accounts, in large part, for the systematic differences we observed between the data and the model fits (Figs. 4⇑–6). Any validity that our model might enjoy thus rests on the possibility that it captures the overall influence of agonist diffusion on airway narrowing dynamics.

In conclusion, we have demonstrated, on the basis of a mechanically based computational model, that the kinetics of diffusive transport of agonist across the airway wall can contribute significantly to the time course of airway narrowing in mice. By inference, this raises the important possibility that alterations in the diffusive barrier presented by the airway wall may play a role in some situations of altered airway responsiveness, especially when significant airway remodeling is involved.

## GRANTS

This study was supported by National Institutes of Health Grants HL-87788 and National Center for Research Resources Centers of Biomedical Research Excellence RR-15557.

## ACKNOWLEDGMENTS

The authors acknowledge the technical assistance of Steve Aimi and Min Wu.

## APPENDIX

Some of the model details below have been presented previously (6). For completeness, we include them here along with the novel additions developed for the present study.

#### Modeling the contracting airway.

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. ASM contraction decreases airway radius by pulling against the parenchymal attachments to the outside of the airway wall. This outward pull comes from the transpulmonary pressure (Ptp), which is determined by lung volume under the assumption of a constant tissue elastance, and by the local distortion of the parenchyma, which is assumed to follow the relationship identified by Lai-Fook (17). The inward recoil of the airway wall is determined by its stiffness, which is assumed to arise from a fraction (1 − *k*) of the airway circumference that expands according to the one-third power of Ptp. The remaining fraction, *k*, of the circumference is assumed to be inextensible, where 0 < *k* < 1. Once activated, the ASM follows the classic Hill force-velocity relationship that is hyperbolic when active force (F_{A}) is less than F_{0}, and linear when F_{A} ≥ F_{0} with slopes matched at F_{0} (13) thus
*r* is airway radius, and *a* and *b* are constants. Following experimental findings reported in rats (9), we set *a* = F_{0}/4. *Equation A1* thus contains two free parameters, F_{0} and *b*.

At any point in time, F_{A} is the active contractile force that adds to the outward recoil of the parenchyma and the inward recoil of the airway wall to give a net force difference of zero. The explicit expression for F_{A} that this produces is derived in Ref. 6, and is given by
*r*_{TLC} is the radius that the virtual hole occupied by the airway would have at total lung capacity (TLC), if it expanded like the rest of the parenchyma, Ptp_{TLC} is Ptp at TLC, and P_{0} is the value of Ptp at which the unconstricted airway induces no distortion in the parenchyma surrounding it. The model fit is insensitive to the value of P_{0}, so we set this parameter in the present study to be equal to the PEEP value of 3 cmH_{2}O used experimentally.

We use the above equations to calculate how *r* varies with time when the model is driven with a prescribed volume signal that, when multiplied by lung elastance, produces a time-varying Ptp(*t*) signal. The volume signal varied sinusoidally with an amplitude of 0.20 ml above the lung volume set by PEEP of 3 cmH_{2}O at a frequency of 200 breath/min. Initially, the ASM is relaxed so that F_{A} = 0, and *r* is determined by the force balance between the inward recoil of the airway wall and the outward recoil of the surrounding parenchyma. Once the ASM in the model is activated, F_{A} is given by *Eq. A2*. The result is substituted into *Eq. A1* to provide d*r*/d*t*, which is then used to determine *r* at the next time step using first-order Euler integration. This new value of *r* is used in *Eq. A2* again to determine the next value for F_{A}, and so on. Finally, invoking the assumption of Poiseuille flow through the airway, a normalized Raw profile is calculated by raising *r*_{TLC}/*r* to the fourth power. This resistance is scaled so as to match the experimentally determined Raw at baseline (precontraction).

#### Modeling methacholine transport.

The transport component of the model begins with an aerosol delivery profile *a*(*t*) as input to the airway wall compartment, where
*T* is the duration of aerosolization. The concentrations of agonist in the airway wall and tissue compartments are governed, respectively, by the following two mass-balance equations:
*k*_{dif} and *k*_{sink} have units of s^{−1}. The values of the parameters of the Hill equation (*Eq. 2*) used in the model at any instant in time are determined by the current value of C_{t} from the nominal values of F_{0}, *a*, and *b* determined from our laboratory's previous study (12) to within an adjustable scale factor γ. That is:

#### Modeling a thickened airway wall.

The output of the model is Raw calculated assuming laminar flow through the airway. That is
*A* is a scaling constant chosen to force the predicted and measured values of Raw to match at baseline (precontraction), and *r* is determined as a function of time as described above. This results in Raw in *Eq. A7* having units of cmH_{2}O/ml.

Let *r*_{0} and *d*_{0} be the radius and airway wall thickness at baseline, respectively, in which case Raw becomes
*r*(*t*) and solving for *d*(*t*), this provides the time course of the luminal radius used to calculate Raw as a function of *t*.

- Copyright © 2012 the American Physiological Society