## Abstract

Recent advances have revealed that during exogenous airway challenge, airway diameters cannot be adequately predicted by their initial diameters. Furthermore, airway diameters can also vary greatly in time on scales shorter than a breath. To better understand these phenomena, we developed a multiscale model that allowed us to simulate aerosol challenge in the airways during ventilation. The model incorporates agonist-receptor binding kinetics to govern the temporal response of airway smooth muscle contraction on individual airway segments, which, together with airway wall mechanics, determines local airway caliber. Global agonist transport and deposition are coupled with pressure-driven flow, linking local airway constrictions with global flow dynamics. During the course of challenge, airway constriction alters the flow pattern, redistributing the agonist to less constricted regions. This results in a negative feedback that may be a protective property of the normal lung. As a consequence, repetitive challenge can cause spatial constriction patterns to evolve in time, resulting in a loss of predictability of airway diameters. Additionally, the model offers new insights into several phenomena including the intra- and interbreath dynamics of airway constriction throughout the tree structure.

- asthma
- intrabreath dynamics
- interbreath dynamics
- multiscale model
- heterogeneity

asthma is a chronic airway disease characterized by reversible hyperreactive bronchoconstriction due to the contraction of airway smooth muscle (ASM) in response to environmental factors. Typical indicators of asthma are increased airway resistance (*R*_{aw}), decreased forced expiratory volume in 1 s, or decreased peak expiratory flow. A standard assessment of lung function in asthmatic patients, the bronchial challenge, is performed by administering an aerosolized ASM agonist and observing the resultant changes in one or more of these indexes. However, recent advances in imaging have revealed that during exogenous airway challenge, lung dysfunction exhibits complex behavior, including a heterogeneous “patchy” distribution of regional ventilation (12, 22, 44) and constriction (7, 49). A long-term clinical study (21) has also shown increased peak expiratory flow variability with loss of temporal correlations in asthmatics, indicating a change in the underlying physiology of asthmatic lungs.

While there have been studies (51, 52) advancing the understanding of the complexities associated with the asthmatic lung, the nature of temporal and spatial heterogeneity of constriction patterns is not well understood even in the normal lung. Indeed, it has been shown that airway diameters after challenge can not be predicted by their initial diameters (9). Furthermore, airway diameters appear to vary greatly on short times scales that cannot be explained by wall remodeling. For example, the pattern of regional ventilation can change from breath to breath, and it is likely that the heterogeneity of flow distribution deep in the lung also changes within a single breathing cycle with unknown consequences for asthmatic attacks.

The spectrum of observed complex behavior may result from the interaction of the increased contractility of the ASM with many apparently unrelated phenomena such as airflow dynamics and agonist distribution as well as pathologies including airway closure, airway inflammation, tissue remodeling, or compromised airway epithelium. Such multiscale problems are difficult to evaluate experimentally, especially in human subjects. On the other hand, current airway models exploring airway behavior are limited in that they do not incorporate dynamics on multiple time and/or spatial scales such as intra- and interbreath dynamics, particle deposition during heterogeneous flow conditions at the level of full airways, or binding kinetics of the agonists at the level of individual ASM cells. Therefore, modeling the multiscale complexity of airway behavior that could be used to better understand the unpredictable nature of asthma attacks still remains unexplored.

In this study, we present a multiscale model that allowed us to simulate aerosol challenge in the human airways during ventilation. We introduce binding kinetics to govern the temporal response of ASM contraction on individual airway segments, which, together with airway wall mechanics and flow conditions, determines local airway caliber. Our model couples local agonist transport and deposition with pressure-driven flow, linking airway challenge with global flow dynamics at time scales ranging from a fraction of a breath to an entire airway challenge.

## METHODS

We constructed a multiscale model to simulate complex spatiotemporal patterns in airway constriction during ASM challenge in normal airways. The airway tree consisted of airways organized into a bifurcating tree that terminated in individual compartments, representing alveolar wall tissue. A neural tone was used to introduce smooth muscle tension, resulting in a slightly heterogeneous distribution of baseline airway diameters. An airborne agonist was delivered to the airways via gas flow, driven by cyclically varying pressure at the tree entrance. In each airway, agonist deposited from the lumen onto a semipermeable epithelial layer and diffused into the wall and surrounding ASM. The time course of individual airway constriction was governed by the kinetics of agonist binding and unbinding to ASM cell surface receptors. Airway caliber was then determined by equations described by Lambert et al. (28) based on morphometric data relating human airway diameters to transmural pressure. Constriction of each airway segment in turn changed its resistance, altered the flow pattern, and hence redistributed agonist delivery throughout the tree. Airway constriction was simulated resulting from a bolus inhalation challenge during pressure-driven ventilation.

### Airway Tree Structure

Our model consists of individual airways organized into a symmetrically bifurcating tree of 10 generations. Using a previously described algorithm (27), unconstricted inner airway diameters (*d*) and lengths (*L = L*_{0}*d*) were assigned for each airway segment, where *L*_{0} is a constant. Wall tissue was treated as incompressible with a fixed volume. Airways maintained a fixed length, and thus outer airway wall diameters (*d*_{o}) were calculated assuming a constant wall area fraction (*f*_{wall}) for all segments as follows:
*A*_{wall} is the wall area and *A* is the area of the lumen.

The tree terminated in 1,024 alveolar compartments, which were modeled as identical spherical sacs with a linear alveolar pressure (P_{a})-alveolar volume (*V*_{a}) relation and alveolar compliance (*C*_{a}) attached to each terminal segment in the last generation of the airway tree:
_{a} is the volume of an alveolus at functional residual capacity (FRC). Gas compression was neglected, and the total alveolar volume at end expiration (P_{a} = 0 for all alveoli) represented FRC and was set to a representative value of 2.5 liters. Total equivalent compliance (*C*_{eq} = Σ*C*_{a}) in normal subjects is ∼0.1 l/cmH_{2}O, tidal volume is 0.5 liters, and airway deadspace volume is 0.15 liters. To meet these constraints, diameters were scaled such that the total airway volume equaled the airway deadspace volume [(π*L*_{0}/4)Σ*d*^{3} = 0.15 liters] and Vfrc_{a}= 2.5 l/1,024.

### Flow Calculations

We assumed laminar flow through the airways and that gas compression and inertia are not significant at the frequency of breathing. Thus, flow (Q = dV/d*t*) through an airway was simply approximated as the ratio of axial pressure drop, the difference in pressures at the proximal end (P_{prox}) and distal end (P_{dist}) of the pipe, to the Poiseuille flow resistance (*R*) of the segment, as follows:
*R* is given by
_{g} is the viscosity of the resident gas at 37°C and 100% humidity. Q is thus generally positive during inspiration and negative during expiration.

Total equivalent resistance of the airway tree (*R*_{aw}) was calculated by adding segment resistances in series and parallel according to the tree structure. Breathing was driven by sinusoidal variation of the external pressure (P_{ext}) at the proximal end of the tree root. When all alveolar compartment pressures are equal (P_{a} ≡ P_{alv}), total flow entering the airway tree can be calculated by substituting P_{ext}, P_{alv}, and *R*_{aw} for P_{prox}, P_{dist}, and *R*, respectively, into *Eq. 3A*. However, heterogeneity causes individual alveolar compartments to possess unique pressures P_{a}, complicating the calculation of individual flows throughout the tree during simulated breathing. Therefore, we introduced the concept of an equivalent pressure (P_{eq}), which acts as the effective pressure at the exit of a segment generated by all downstream airways and alveolar compartments. As detailed in the appendix, this formulation allows the treatment of nonhomogeneous alveolar pressures. Starting at the terminal ends, we recursively calculated P_{eq} to obtain P_{eq} of the root of the tree (P_{eq,aw}) and to obtain the total flow into the tree (Q_{tot}) from *Eq. 3A* using P_{ext}, P_{eq,aw}, and *R*_{aw}. Next, we substituted P_{eq,aw}, flow into the tree (Q_{in}) = Q_{tot}, and *R* of the root back into *Eq. 3A* to determine the actual P_{dist} of the root. Realizing that P_{dist} of the root equals the proximal pressure of both daughter segments of the root, we can sequentially resolve the instantaneous flow through all segments by repeated application of the P_{eq} function for each airway until we reach the alveolus, where, for the last generation segment, P_{dist} = P_{eq} = P_{a}. P_{a} was then updated as follow:
_{a}) at *generation 10* is the exit flow from the terminal airway, which is attached to the alveolar compartment.

### Agonist Transport and Deposition

During simulated aerosol challenge, air entering the airway tree contained a prescribed concentration of airborne ASM agonist particles that traveled through each airway with the local flow (see Fig. 1*A*). Particle size was assumed to be small, in the range where gravitational and inertial effects on the particles are negligible (18). Furthermore, to couple flow with agonist deposition in a large tree, we assumed that the agonist concentration (*X*) was homogenously distributed within each airway and thus neglected intra-airway axial and radial concentration gradients. Generally, the amount of agonist (d*S*) passing through a cross section of a segment parallel to the flow (Q) in a short time (d*t*) is simply the product of Q through that plane and the concentration of agonist in the air:
_{in} is the flow entering a given airway segment (either from the proximal end during inspiration or the distal end during expiration) and *X*_{in} is the concentration of agonist at the bifurcation from which the gas entered (the parent airway for inspiration and the daughter airways during expiration), then the rate *S*′_{in} describes the amount of agonist entering the airway per unit time, as given by
_{out} = Q_{in} at a rate *S*′_{out}, as given by
*X*_{aw} is the agonist concentration in the airway. A fraction (*S*_{dep}^{′}) of the agonist in the airway lumen was assumed to deposit onto the epithelial layer per unit time, a semipermeable layer of cells separating the lumen from airway wall compartment, proportional to the agonist concentration in the airway (*X*_{aw}) and the surface area of the airway lumen (*A*_{s}), as follows:
*S*′_{aw} is given by the following:
*X*_{bif}) was computed as the flow-weighted average of the incoming agonist concentrations as follows:
*X*_{bif} simplifies to *X*_{aw} of the parent airway and was taken as the *X*_{in} for the two daughters. During expiration, *X*_{bif} was the weighted sum of *X*_{aw} values corresponding to the two daughters and hence was the input *X*_{in} for the parent. We note that *Eq. 8* is also valid during pendelluft effects. Finally, the amount of agonist that reaches the alveoli by the end of inspiration does not deposit onto the alveolar walls; instead, the alveolar concentration of the agonist (*X*_{a}) was taken as the input *X*_{in} for the corresponding terminal airways during the next expiration (see the schematic representation shown in Fig. 1*A*).

### Modeling ASM Challenge

#### Agonist diffusion and binding kinetics.

We developed a simple pharmacokinetic model to describe the time course of agonist administration in a single airway segment. The model consists of a single compartment representing both the airway wall tissue and ASM (Fig. 1*B*). Within the tissue compartment, the agonist can exist in one of two states: either bound (B) or unbound (U) to the ASM cell receptors, with amounts *S*_{B} and *S*_{U}, respectively. Free agonist molecules entered the tissue compartment through deposition onto and diffusion through the epithelial layer surrounding the airway lumen at a rate *S*_{dep}^{′}. Unbound agonist binds to ASM receptors at a rate α, whereas bound agonist unbind at rate β. Unbound agonist molecules were also continually removed from the tissue compartment by diffusion into the pulmonary circulation at a rate γ. Assuming a steady flow of blood, agonist entering the vasculature was unable to return to the airway tissue. Receptor-bound agonist molecules were unable to directly enter the bloodstream, but rather had to first unbind with rate β to reenter the unbound state. The differential equations governing the rate of change of *S*_{B} and *S*_{U} are as follows:
*t*βγ and χ_{B} = (βγ/α)*S*_{B}. This equation is exactly that of a damped oscillator with natural frequency (ω_{n}) = βγ and damping coefficient (ζ) = (α + β + γ)/2βγ. Since the physiologically relevant range of our parameters is α, β, γ > 0, the system above corresponds to an overdamped (ζ > 1) oscillator. Changes in receptor binding kinetics can be investigated by varying α, β, or γ.

#### Neural input.

To account for vagal stimulation of the ASM, a small amount of agonist was locally applied to each airway independent of the inhaled agonist. This neural agonist mimics the cholinergic input from the parasympathetic nerves that innervate the ASM but is assumed to act in the same manner as our inhaled agonist. The level of neural stimulation (*S*_{N}) was assigned to each airway from a uniform random distribution, allowing us to introduce a slight heterogeneity in airway caliber into an otherwise perfectly symmetric tree.

#### Force response of the ASM.

Experimental studies (36, 53) of ASM strips have described a sigmoidal log-dose response curve to several agonists, which is characteristic to many biological response phenomena (15). In these studies, ASM cells were placed in fluid baths of varying agonist concentrations until equilibrium contractile force was achieved, and the force generated was plotted as a sigmoidal function of the logarithm of the agonist concentration. However, in general, the agonist concentration can be different in the airway lumen than in the airway tissue, and the difference might vary along the breathing cycle as well as from cycle to cycle. Thus, in terms of our pharmacokinetic model, it is the actual number of bound agonist-ASM receptor complexes that determines force generation. Accordingly, the tension generated by the ASM (T_{ASM}) can be described as a function of the instantaneous number of agonist bound receptors as follows:
_{ASM} is the reactivity (maximum achievable tension) of the ASM and B_{ASM} is the sensitivity (dose required to achieve 50% of maximal tension)

### Airway Segments

We used Laplace's law (pressure = 2T_{ASM}/*d*), assuming negligible wall thickness, to convert ASM tension on the outer airway wall to a pressure and then subtracted this from the gas pressure (mean of the proximal and distal airway pressures) to obtain an effective transmural pressure. Airway caliber was calculated using the equations proposed by Lambert et al. (28), which describe the equilibrium airway diameter as a function of the effective transmural pressure.

### Simulations

#### Single airway simulations.

We simulated a bolus challenge of a single isolated tracheal airway for four parameter selections of our two-state binding kinetics model. With the analytic solutions of the impulse response of the two-state model (Fig. 2*A*), and our airway pressure-diameter relation described above (Fig. 2*B*), we were able to predict the change in resistance of a single airway due to an impulse of agonist during simulated breathing (Fig. 2*C*).

#### Airway tree simulations.

All further simulations were performed for a 10-generation tree composed of 2,047 airways and 1,024 alveolar compartments. The total equivalent resistance of the tree was used as a measure of *R*_{aw}.

#### Baseline R_{aw} variability.

To study the baseline (including tone but excluding challenge) variability of *R*_{aw} from neural ASM tone, the tree was first initialized by randomly selecting *S*_{N} values for each airway, calculating the resulting equilibrium diameters, and then reevaluating *R*_{aw} normalized by its initial (pretone) value. This was repeated 1,000 times while different values for *S*_{N} were drawn from a specified uniform random distribution. Kernel smoothing (Gaussian) was applied to the histogram of *R*_{aw} values to estimate its probability density distribution.

#### Inhaled agonist challenge.

Baseline neural tone was implemented once as described above and fixed throughout the remaining simulations. Breathing was simulated by applying a 0.25-Hz sinusoidal pressure waveform [P_{ent}(*t*); amplitude = 2.5 cmH_{2}O] to the tree entrance for the duration of 50 breaths, resulting in a 0.5-liter tidal volume and 7.5 l/min ventilation prechallenge. To simulate an airway challenge maneuver, a constant input agonist concentration was inhaled during the inhalation phase of the second breath. The system was numerically solved with a 0.001-s step size using MATLAB. Airway segment diameters (*d*) and *R*_{aw} were recorded as a function of time every 10th solution step. Table 1 shows the parameter values of the model.

#### Data analysis.

The above variables were considered both at baseline and at the time instant when *R*_{aw} reached its peak value. Additionally, pathway resistances (ρ_{a}), defined as ρ_{a} = (P_{in} − P_{a})/Q_{a}, were calculated directly from airway segment resistances at the corresponding time points for all 1,024 terminal segments (see the appendix). We investigated three specific potential sources of airway hyperresponsiveness: a compromised epithelial layer (δ), an increased ASM agonist binding rate (α), and an increased unbinding rate (β). Specifically, we determined the sensitivity of *R*_{aw} and constriction heterogeneity (taken as the SD of pathway resistances) to variations in model parameters.

## RESULTS

Figure 2 shows the pharmacokinetics of a single airway segment representing the root of the airway tree in response to an impulse challenge. Figure 2*A* shows the analytical solutions to the two-state model described by *Eq. 9* for *S*_{B} in time for various parameter values (see Table 1). Doubling α, which corresponds to increasing the binding rate of an agonist to a receptor, amplified the number of bound receptors for all times, with a minimal effect on the time to reach peak constriction. Doubling β, which corresponds to increasing the unbinding rate of an agonist to a receptor, decreased the number of bound receptors while also resulting in a shorter time to reach a peak. Scaling all parameters by a constant factor decreased the rise and fall times while not influencing the peak number of bound receptors. Note that due to the symmetry in β and γ (see text after *Eq. 10*), varying γ is identical to varying β and thus was omitted from the analysis. The tension generated by the ASM layer was calculated from the number of bound receptors using *Eq. 11*. Figure 2*B* shows the instantaneous pressure-area relation of the “intial” airway segment from Fig. 2*A* at different time points before challenge (*curve a*), at half peak (*curve b*), and at the peak number of bound receptors (*curve c*). The lumen area approached zero as pressure decreased to negative infinity and saturated at high positive pressures. As the number of bound receptors changed during a simulated challenge, the pressure-area relation constantly shifted in time, as shown in Fig. 2*B*. Figure 2*C* shows the combined effects of binding kinetics and pressure-area characteristics on the resistance of the segment during a simulated agonist challenge. The number of bound receptors followed the results shown in Fig. 2*A* and the lumen area followed the curves shown in Fig. 2*B* while the segment was driven by a sinusoidal pressure input. Instantaneous resistance of the segment was then calculated directly from the lumen area.

To explore the effect of neural tone on *R*_{aw}, each of the 2,047 airways were randomly assigned a tonic input *S*_{N} from a prescribed uniform distribution. The different realizations of the random input produced a distribution of possible *R*_{aw} values. Figure 3 shows the estimated *R*_{aw} probability density distributions for different ranges of the tonic input. For the smallest tone intensity, the *R*_{aw} probability density distribution was narrow, with a mean of ∼120% of the resistance of the tree pretone. Increasing the neural tone amplified the mean and widened the *R*_{aw} distribution. The smallest range (0.1–0.2) was used for initializing airway trees for all further simulations.

The configuration of a section (*generations 4–10*) of the airway tree at peak constriction, with the parameters shown in Table 1, is shown in Fig. 4*A*. From the imposed neural tone and subsequent interference of the flow pattern, the agonist was unevenly distributed throughout the tree and generated a heterogeneous distribution of segment diameters. A second inhalation challenge, after recovery from the first, produced nearly identical peak *R*_{aw} despite a different distribution of diameters (Fig. 4*B*), implying the existence of memory in the system. Notice that many regions were either more or less constricted than in the first simulation, despite the same neural tone and agonist dose (see arrow).

Figure 5 shows the change in diameters over time in four airways randomly selected from *generation 6* of the tree. Each trace was normalized by the baseline diameter after the application of neural tone. The small asymmetry due to neural tone accounts for large variations in the behavior of airways of the same generation. Note that the least constricted (dashed lines) airway at peak challenge (time ∼ 35) remained constricted for a longer duration than a more constricted airway (solid line). However, an intermediate airway (dotted line) constricted nearly as much as the latter but remained constricted slightly longer than some less-reactive airways. One airway (dashed-dotted line) constricted over 50% further than the others, despite nearly identical initial conditions. All airways eventually reached equilibrium at their baseline diameters, as the inhaled agonist was entirely washed out of the ASM compartment.

The diameters of all airways at peak constriction (normalized by their baseline diameters) were related to their baseline diameters (Fig. 6*A*). Each point corresponds to 1 of the 2,047 total airway segments. The patterns of bands were formed because airways at each generation had almost identical initial diameters with the small variability caused by the random neural tone. Notice that this slight asymmetry across each generation caused a substantially greater variation in diameter after peak constriction was reached, making the prediction of constricted diameter from baseline conditions difficult.

The 1,024 individual pathway resistances after peak constriction was reached were plotted against their baseline values before challenge (Fig. 6*B*). Pathways that were initially of low resistance tended to increase to a greater extent than those with an initially high resistance. When the results shown in Fig. 6, *A* and *B*, were compared, it appeared that the correlation between pre- and postconstriction was weaker for individual diameters than pathway resistances. Although segments that lie on the direct path from the root to the specific terminal airway contributed the most to its pathway resistance, each pathway resistance still depended on the resistance of every other airway in the tree. Pathways with low resistance received increased flows and higher rates of deposition. Thus, the segment's diameter after peak constriction was reached correlated better with the baseline pathway resistance than its own baseline diameter, which we attribute to the coupling of agonist delivery with local flows.

The sensitivity of the model output to several parameters is shown in Fig. 7. The sensitivity curves correspond to the peak constriction and are normalized so that unity corresponds to the baseline parameters. *R*_{aw} monotonically increased with increasing binding rate and the rate of absorption of the agonist, whereas it decreased with increasing unbinding rate. The heterogeneity of the constriction pattern characterized by the SD of the pathway resistances showed different patterns for each parameter. The decline of SD when a parameter was further increased is due to the fact that once about half of the airways were constricted; further constriction increases *R*_{aw} but reduces heterogeneity.

## DISCUSSION

In this study, we introduce a model that simulates dynamic pressure-driven breathing in an airway tree that incorporates tissue properties, flow dynamics, and agonist distribution together with pharmacokinetics over a wide range of time scales. To our knowledge, no similar model has been proposed. The integration of these specific components allowed us to better comprehend how small-scale phenomena, such as binding kinetics, impact large-scale dynamics, such as the flow pattern over the entire tree, producing observable global functional characteristics (e.g., *R*_{aw} and constriction heterogeneity). The model also captures intrabreath dynamics while producing a time series of hundreds of breaths, linking pressure and flow fluctuations within a single breath to the entire time course of an inhalation challenge. In addition to presenting this model, we also introduced two new concepts: P_{eq} and pathway resistance, which not only allow for an efficient calculation of intrabreath dynamics but also provide a better intuitive understanding of constriction dynamics in a bifurcating tree structure.

### Assumptions and Limitations

Breathing involves complex interactions of tissue mechanical properties, fluid flow, particle deposition, gas exchange, and the cellular response, each of which are computationally expensive to simulate and, when taken together, prove to be a significant computational challenge. Thus, it was necessary to invoke several carefully chosen assumptions to produce a computationally practical model while preserving the ability to simulate dynamic breathing over several length and time scales.

Viscous fluid flow through a bifurcating tree structure can be highly nonlinear and depends on the local Reynold's number, boundary conditions, branching angle, oscillation frequency, and flow profile distortions (40). A further complication arises from the fact that ASM challenge involves particle deposition, which is influenced by inertial, gravitational, and diffusive effects depending on the particle shape and size distribution (17, 23). Solving for the flow and deposition alone through a tree of even a few generations is computationally intense, usually requiring finite-element or computational fluid dynamic techniques (37, 54), and is completely impractical for a 10-generation tree. While the implementation of asymmetric bifurcations is possible, adopting a symmetric one-dimensional tree allowed us to distinguish between the effects of structural asymmetry and flow dynamics on constriction heterogeneity. Thus, we constructed a symmetric airway tree in which each airway segment had only one degree of freedom (radius) and thus lacks spatial information such as branching angle, axial and circumferential heterogeneity, and positional dependence on gravity.

At a lower level, the binding kinetics were assumed to obey linear dynamics, where we described binding, unbinding and washout with simple rate constants that were independent of agonist solubility in the tissue and blood and changes in bronchial perfusion. We also discounted other potential rate-limiting steps in the ASM force response, such as the finite number of receptors and intracellular second messenger systems. The rates are difficult to asses experimentally since each rate constant describes a combination of events such as diffusion into the ASM lipid membrane, repetitive binding of active sites of agonists attached to secondary sites on the receptor, enzyme-driven degradation, diffusion into the bloodstream, etc. It would require a systematic and highly involved series of experiments to separate out each individual rate constant in vivo, and this is outside the scope of the present study. Accordingly, our parameters were selected to qualitatively agree with the time course of airway constriction seen in both humans and animals rather than model any specific agonist and tissues properties.

Airway caliber has several sources of time dependence, which are generally difficult to separate. For instance, dynamic airway constriction results from ASM force generation (35) and lumen pressure fluctuations (43) as well as the viscoelastic nature of the ASM and surrounding parenchymal tissues (34). In the model, pressure and flow fluctuations arise from the input driving pressure, which is at a fixed frequency. The ASM response is on a slightly longer time scale of 10 breaths, depending on the parameter values of the binding kinetics. Tissue viscoelasticity is known to follow a power law with even longer relaxation times (6, 47), making it difficult to distinguish between changes in airway diameter resulting from pressure or the agonist. For this reason, and also to further reduce computational complexity, we assumed that mechanical equilibrium in the tissues was instantaneous. Thus, the time scale of pressure-induced fluctuation in our model was well separated from that of ASM force generation. Although this may limit the accuracy of our model results on short time scales, we were primarily interested in the dynamics during the entire course of an airway challenge. We note, however, that it is possible that tissue viscoelasticity leads to airways reaching a dynamic equilibrium (30), which is not necessarily the static equilibrium determined from our model (3, 4). Additionally, passive wall or smooth muscle properties may also change with pressure fluctuations and ASM activation (29, 45). Furthermore, we also neglected the direct effects of tidal stretch on ASM contractility (16, 19) and load dependence of shortening velocity. These remain limitations of the present model. Nevertheless, this simpler network still allowed us to explore the dynamics of a system at a size (2,047 segments) otherwise impossible.

### Agonist Transport and Flow Coupling

Airway smooth muscle contraction begins with the binding of specific molecules to their corresponding cell surface receptors. This induces an intracellular signaling cascade, resulting in active force generation within the cell (46). While the corresponding deformation is limited to the surrounding airway tissue, local changes in airway caliber still alter the global flow pattern because segment constriction diverts air flow to other regions of the lung. The agonist is transported throughout the airway tree along with this flow; thus, local agonist-induced constriction is coupled to agonist transport. This mechanism establishes a negative feedback between agonist deposition, ASM constriction, and gas flow, which, under normal conditions, would prevent highly constricted airways from receiving further agonist. This coupling has two specific implications. Under normal circumstances, modest airway constriction acts to distribute inhaled agonists evenly throughout the airway tree. This would essentially reduce, rather than amplify, small heterogeneities in airway caliber in the presence of irritating airborne agonists. In terms of airway function, a broad and relatively homogeneous constriction of all airways generates lower total *R*_{aw} and ventilation heterogeneity compared with regionally clustered severe constrictions (33). Indeed, regions that are initially exposed to slightly more agonist than the average airway may overreact, becoming severely constricted and diverting more agonist into neighboring airways. The neighboring regions would then receive higher exposure to the agonist and constrict, while the previously constricted regions would be able to relax. We speculate that this protective effect may be strongly present in normal subjects. In asthma, however, such a stabilizing effect of agonist/flow coupling is likely to be disrupted through ASM hyperreactivity or some other means, such as the instabilities proposed by Anafi and Wilson (2) and more recently by Affonce and Lutchen (1). As a consequence, clustered constriction patterns can suddenly form due to the presence of even the smallest heterogeneity in the tree (51), a feature reminiscent of chaotic systems. If this is the case, a second implication of flow-agonist coupling is that clusters of severely constricted airways can evolve spatially in time in the presence of long-term agonist exposure.

### Constriction Dynamics

As shown in Fig. 5, we demonstrate that airways of a single generation can exhibit different reactions during an aerosol challenge. In particular, some airways constricted to a greater extent than others. It also appears that some airways reached baseline diameters after challenge quicker than others. This is due to variations in neural tone and agonist deposition, along with the sigmoidal dose response of ASM force generation. An ASM cell that has already reached its upper plateau of force generation cannot constrict much further. However, there would be a larger apparent reaction if that cell was in the linear region of its dose-response curve and received the same amount of agonist. On the other hand, the length of time the airway remains constricted depends on the total deposited amount of agonist. Thus, initially unconstricted airways may remain constricted longer than preconstricted airways. Upon subsequent challenge, the initial neural tone distribution will, in general, be different, and the global flow pattern will deviate from the original challenge. The fact that prechallenge diameters vary due to differences in structure, baseline neural tone, or prior agonist exposure makes it nearly impossible to correlate postchallenge diameters to prechallenge diameters without knowledge of the complete history of airway diameters and ASM force generation as well as agonist deposition.

In this study, we simulated short-term, concentrated airway challenges, similar to bronchial challenge maneuvers performed in normal subjects. Our results indicate that challenge of an almost symmetric airway tree can result in topologically correlated constriction patterns. However, we did not find patchy constriction and ventilation patterns (Fig. 4*A*), as suggested by other models (51, 52) as well as by PET (22) and MRI (50). This is likely due to the above-mentioned protective effect of the negative feedback between agonist transport and constriction. The only heterogeneity in our model was introduced through the neural tone, which alone does not lead to patchy ventilation. However, there may be other sources of baseline heterogeneity, such as smooth muscle thickness or epithelial damage, that could work against the protective negative feedback of agonist flow coupling. It is possible that combining intrabreath dynamics with airway instabilities that lead to airway closure would lead to such gross heterogeneity of ventilation.

### Two-State Pharmacokinetics

The time course of local airway constriction during challenge is dependent on the kinetics of ASM contraction, which are fundamentally determined by molecular-level interactions (10). We therefore implemented a simple kinetic model of ligand-receptor binding on the ASM surface. Our ligand, the ASM agonist, can bind to ASM receptors to form a receptor-ligand complex with a certain probability that determines the rate constant. The formation of bound complexes induces ASM force generation, such as methacholine binding to the muscarinic M_{3} receptor (10). The complex dissociates at a different rate, which depends on the binding affinity of the ligand to its receptor. Furthermore, ligand molecules can be removed from the vicinity of the receptors through diffusion, degradation, or other processes at a third lumped rate. This system, although very simple, can account for lumped molecular processes. When built into a tree structure, the combined model can provide an organ-level transient response to challenge with qualitative agreement with physiological data (5, 8). In some respect, our model may appear to be similar to that used by Lauzon and Bates (31); however, they only considered the transport of agonists but not the separate binding and unbinding processes.

### Neural Tone

It is well established that smooth muscle tissue is innervated by parasympathetic nerves (13, 41). These nerves can prompt cholinergic ASM contraction when activated (14). In animals, it has been observed that this neural regulation plays a primary role in establishing a baseline ASM tone, which was shown by the marked reduction in airway resistance when the vagal nerve was severed (11, 42). It is feasible to assume that this phenomenon also exists in humans (32). Therefore, we attempted to incorporate neural regulation in our baseline airway model. As specific data on neural regulation at the level of individual airways are not well known, we assumed that neural activation varies randomly within a small range. Since the regulation is cholinergic in nature, we reasoned that neural activation acts similarly to our agonist through a receptor-ligand binding process (13). It is likely that neural regulation does change over time, so the same airway tree would achieve different realizations of neural activation measured at different occasions. In our simulations, the different possible realizations of neural activation produced a distribution of total *R*_{aw} (Fig. 3), and this distribution may correspond to the intrasubject variation of *R*_{aw}. Changes in neural activation, and thus baseline ASM tone, may also account for the poor temporal correlation of airway diameters found recently in dog lungs (25). The effect of challenge, in which flow and agonist distribution is coupled and superimposed on the neural tone, is to break the correlation between initial diameter and constricted diameter (Fig. 6*A*), in good agreement with recent computed tomography imaging results in dogs (9).

### Pathway Resistance

We defined pathway resistance as the ratio of the flow through a terminal airway segment to pressure drop from the root to the end of that segment. Experimentally, pathway resistance can be measured using the alveolar capsule technique (20), and it relates precisely to the time constants of inflation of the alveolar compartments terminating each path. Thus, the distribution of pathway resistance correlates with the distribution of alveolar tidal volumes, providing the link between heterogeneity in airway constriction to the resulting ventilation profile. On the other hand, heterogeneity in alveolar compliance while holding pathway resistances equal has a similar effect of imposing heterogeneous ventilation (38). As a consequence, maintaining an even distribution of ventilation throughout the lung requires an inverse relationship between pathway resistance and alveolar compliance and, hence, a narrow distribution of pathway time constants. In the normal lung, this should be related to the design of an optimal structure for gas exchange (27). Indeed, acinar volumes show significant variation (39) similar in magnitude to the variability of alveolar pressures using the capsule technique (20). In diseases, however, the time constant distribution is significantly broadened (33).

During heterogeneous airway constriction, it remains difficult to make accurate predictions of changes in airway and tissue properties from measurements at the airway opening. Consequently, distributed airway models have been developed and applied to observations of lung input impedance in the presence of airway challenge to partition lung tissue responses from inhomogeneous airway constriction (25, 26, 48). These models include a set of parallel pathways with a continuous distribution of resistances (48) or viscoelastic tissues (24). Indeed, another way to interpret pathway resistances is via the graph transformation of a tree with *N*_{t} terminal edges into a parallel network of *N*_{t} resistors. The resistances in the parallel network correspond directly to the *N*_{t} pathway resistances of the tree. Since we have shown that pathway resistances are not independent, the parallel resistances are also not independent. Nevertheless, our calculations allowed us to determine the actual distribution of *N*_{t} parallel resistances from any arbitrary tree structure and may enhance the utility of the distributed models. Specifically, we found that a uniform heterogeneity in airway tree diameters (due to neural tone) created a slightly skewed, unimodal pathway resistance distribution. However, constriction heterogeneities higher up in the airway tree greatly distort this into a multimodal distribution. This implies that the pathway resistance distribution may also be useful in differentiating a diffuse constriction present throughout the airway tree from a serial constriction of individual airways of relatively large diameter.

### Conclusions

In summary, we introduced a novel model that combines macroscale flows and agonist transport with microscale binding kinetics to predict ASM contraction in a dynamic manner. The model reveals a novel protective characteristic of the normal airway tree and offers new insights into several phenomena, including intra- and interbreath dynamics of airway constriction throughout the tree structure.

## GRANTS

This work was partially supported by National Heart, Lung, and Blood Institute Grants HL-059215 and HL-090757.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

- Copyright © 2010 the American Physiological Society

## APPENDIX

#### Equivalent Pressure (P_{eq})

P_{eq} is a resistance-weighted pressure. To introduce P_{eq}, we started with a two-generation tree, a single bifurcation, and labeled the parent and the left and right daughters as segments (0), (1,0), and (1,1), respectively. For the case of homogenous terminal pressures P_{1,0} = P_{1,1} = P_{out}, the total entrance flow was calculated as follows:
*R*_{eq} = *R*_{0} + (*R*_{1,0}*R*_{1,1})/(*R*_{1,0} + *R*_{1,1}) is the equivalent tree resistance of the network and P_{ent} is the pressure at the entrance of segment (0). The situation is more complex when P_{1,0} ≠ P_{1,1}. However, it is possible to determine the effective homogenous pressure P_{eq} that plays the role of P_{out} in *Eq. A1*. By denoting the exit pressures of segments (0), (1,0), and (1,1) as P_{0}, P_{1,0}, and P_{1,1}, respectively, we solved the following set of four equations with four unknowns for the total flow:
_{0} = Q_{ent}. We then arrived at the following solution for a two-generation tree:
_{eq} is the effective pressure at the exit of segments (1,0) and (1,1) that can be used together with the total *R*_{eq} of the bifurcation to obtain the flow entering segment (0). Notice that the calculation of P_{eq} does not require the knowledge of flows, a feature that allows for an easy calculation of individual flows in a large tree when terminal pressures are heterogeneous. Consider, for example, a three-generation tree. The computations start at the bottom with four different terminal pressures at *generation 2*. The first step is to obtain two separate P_{eq} values at *generation 1* using *Eq. A4*. These are now taken as the terminal pressures P_{1,0} and P_{1,1}, and, using *Eq. A4* again, P_{eq} of the root can be obtained and the flows calculated as for a two-generation tree.

#### Pathway Resistance

We define the pathway resistance as the ratio of the pressure drop across the path from the tree root to the exit of a given segment to the flow through that segment. For example, the pathway resistance to segment (1,0) was calculated as follows:
_{1,0} that depended solely on segment resistances. From *Eq. A3*, we obtained the same value for Q_{ent} when P_{1,0} and P_{1,1} were replaced by P_{eq} and thus the same P_{0} since P_{ent} and *R*_{0} are constants. Therefore, we could simplify our tree structure to two resistive elements in series: the root with resistance *R*_{0} and the parallel sum of *generation 1* resistances (*R*_{∥}). Now, we used the potential divider relation to simplify *Eq. A5* as follows:
*Eq. A6* into *Eq. A5*, we obtained the following final expression for pathway resistance:
_{1,1} was calculated similarly, by replacing *R*_{1,0} with *R*_{1,1}. The algorithm can be extended to larger trees by successively replacing the terminal segments (e.g., *R*_{1,0}) with new bifurcations.