Abstract
Suki, Béla, Huichin Yuan, Qin Zhang, and Kenneth R. Lutchen. Partitioning of lung tissue response and inhomogeneous airway constriction at the airway opening. J. Appl. Physiol. 82(4): 1349–1359, 1997.—During a bronchial challenge, much of the observed response of lung tissues is an artifactual consequence of inhomogeneous airway constriction. Inhomogeneities, in the sense of time constant inequalities, are an inherently linear phenomenon. Conversely, if lung tissues respond to a bronchoagonist, they become more nonlinear. On the basis of these distinct responses, we present an approach to separate real tissue changes from airway inhomogeneities. We developed a lung model that includes airway inhomogeneities in the form of a continuous distribution of airway resistances and nonlinear viscoelastic tissues. Because time domain data are dominated by nonlinearities, whereas frequency domain data are most sensitive to inhomogeneities, we apply a combined timefrequency domain identification scheme. This model was tested with simulated data from a morphometrically based airway model mimicking gross peripheral airway inhomogeneities and shown capable of recovering all tissue parameters to within 15% error. Application to our previously measured data suggests that in dogs during histamine infusion 1) the distribution of airway resistances increases widely and2) lung tissues do respond but less so than previously reported. This approach, then, is unique in its ability to differentiate between airway and tissue responses to an agonist from a single broadband measurement made at the airway opening.
 lung resistance
 lung elastance
 tissue resistance nonlinearity
 distribution of resistances
 Wiener structure
 Horsfield model
the importance of the contribution of lung parenchymal tissue resistance (Rti) to total lung resistance (Rl) at physiological frequencies (0.1–1 Hz) has recently become evident (15, 9,1113, 1922, 27). Several studies have reported that in response to bronchochallenge Rti and lung tissue elastance (Eti) are significantly elevated (1, 2, 7, 9, 12, 13, 16, 18, 19, 21, 22, 28, 31) and show an increased nonlinear inverse dependence on tidal volume (Vt) delivered (1, 4, 22, 26,30, 3234). Simple partitioning of Rl to tissue and airway components has been attempted in control as well as in challenged conditions (1, 5, 1113, 16, 1923, 27). Bates et al. (2), however, suggested that in dogs, inhomogeneous airway constriction contributes significantly more to the increases in dynamic lung elastance (El) than lung tissues. Therefore, the partitioning of Rl to tissue and airway components can be subject to large errors that may then lead to false conclusions about the responsiveness of parenchymal tissues. Indeed, the influence of inhomogeneities on the partitioning has been confirmed experimentally (21) and inferred via modeling analysis in two other recent studies (16, 20). Nevertheless, to our knowledge there is no method available to reliably partition airway and tissue responses from pressureflow data measured noninvasively at the airway opening in the presence of inhomogeneous airway constriction and nonlinear tissue viscoelasticity.
In a recent study, we addressed one aspect of the above issue. Identifying tissue nonlinearities from data taken at the airway opening, we examined the correlations between the changes in the frequency and amplitude dependences of lung tissue by using a nonlinear blockstructured modeling approach (33). We found that the best lung structure was a linear airway compartment consisting of an airway resistance (Raw) and inertance (Iaw) in series with a nonlinear model for the lung tissues. The nonlinear tissue model was a socalled “Wiener structure,” which is a cascade connection of a linear viscoelastic block and a nonlinear zeromemory block (Fig.1). This model provided excellent fits to the pressureflow data in both control and in the late response to histamine that occurred 15–20 min after the histamine infusion was started. Moreover, our analysis showed that the primary cause for the increased Vt dependence during constriction is not a change in the nonlinear mechanisms. Rather, it results from an exacerbation of the preexisting nonlinear mechanisms through an increase in the magnitude of the linear tissue impedance. This conclusion is also in agreement with some recent histological findings (17).
In our previous study, we fit data during the late response to histamine challenge (33). In the late response, the distribution of alveolar pressures (assessed by the alveolar capsule technique) suggested that airway inhomogeneities had decreased compared with the peak response that occurred in 1–3 min after infusion (22). Nevertheless, we cannot rule out that inhomogeneities also contributed to the estimates of Rti and Eti. Moreover, during the peak response the capsules displayed significant inhomogeneities. This paper addresses the following question: Is it possible to distinguish from pressureflow measurements made at the airway opening whether the lung tissues underwent real changes, or have Rl and El increased only due to the presence of inhomogeneities? From the inputoutput point of view, parallel inhomogeneities are an inherently linear phenomenon, whereas real tissue changes have been found to be accompanied by increased Vt dependence (33, 34). Thus we hypothesized that new model classes that include different types of airway inhomogeneities as well as tissue nonlinearities will allow us to differentiate between airway and tissue responses. We tested our hypotheses by developing and fitting such models to pressureflow data measured previously in histaminechallenged dogs and simulated via an anatomically based structural model.
MODEL DEVELOPMENT
Nonlinear homogeneous model.
First, we summarize our previous nonlinear model, which has been detailed (33). Briefly, the airway tree was modeled by a linear block (Law), whereas lung tissues were represented by a linear tissue block (Lti) and a nonlinear tissue block (Nti). Assuming a homogeneous airway system, the airway compartment was represented by an equivalent tube, having a linear flow resistance Raw and gas inertance Iaw. Thus the pressureflow relationship of the Law block is defined as
The lung tissue is described by a linear viscoelastic subsystem (Lti) and a nonlinear zeromemory subsystem (Nti) arranged in cascade (33). When the Lti precedes Nti, the structure is called the Wiener model (see Fig. 1). We previously showed that during tidal breathing, the Wiener model provides superior descriptions of lung tissues compared with other model structures (33). The model of Lti is motivated by the power law dependence of stress relaxation in a rubber balloon, described by Hildebrandt (14) and adapted by Hantos and coworkers (11,12) for lung tissues. This model describes the linear impedance of the tissues Z_{Lti} as
For the nonlinear tissue block Nti we found that the following thirdorder polynomial provided appropriate balance between quality of fit and number of free parameters
Nonlinear inhomogeneous models.
To modify the above model to account for the presence of airway inhomogeneities (25), we represented the airway tree by a set of airway pathways arranged in parallel (Fig. 2) (see Ref. 12). Each branch is modeled as an inseries connection of an airway resistance, inertance, and a linear tissue block the impedance of which is given in Eq. 3 . To add nonlinearity, we connected this structure in cascade with the nonlinear tissue block described by Eq. 4 (see Fig. 2).
The variation in airway resistance from path to path (Fig. 2) is implemented in a probabilistic manner. If we distribute the pathway resistances according to a distribution functionn(R_{i}), we can calculate the total admittance of the linear part of the system (Y
_{L}), which is the summation of the conductance of the individual branches
By choosing different values for μ, various model classes can be obtained. Because the model is intended to be used in the inverse sense, i.e., fitting data, we are interested in μ values that lead to a closed form solution forY
_{L}. Simple analytical solutions are obtained when μ takes integer values of −1, 0, and 1. When μ = −1,n(R) increases linearly with R; when μ = 0, R is uniformly distributed over all the branches; and when μ = 1,n(R) decreases hyperbolically with R. Finally, the total impedance of the linear system Z_{L} can be calculated by taking the inverse of the admittanceY
_{L} given inEq. 7
. Accordingly, Z_{L} for the three different μ values are obtained as following
METHODS
Description of experiments.
The data used in this modeling study were published earlier (22) and were used in our recent modelfitting exercise (33). Briefly, the data were acquired with an optimal ventilator waveform (OVW), which is a broadband waveform containing energy from 0.23 to 5 Hz (23). The frequencies in the OVW are selected according to a nonsum nondifference (NSND) criterion of order 3, which eliminates harmonic distortion and minimizes crosstalk at the frequencies present in the input flow waveform and thereby provides smooth estimates of the input impedance of the system being measured (32). The spectrum and the phase angles of the components are optimized to ventilate the subject. The OVW approach with a Vt of 300 ml was applied in four anesthetized, tracheotomized, paralyzed, and thoracotomized dogs before and during histamine infusion (16 μg ⋅ kg^{−1} ⋅ min^{−1}). Airwayopening pressure (Pao) with reference to atmospheric pressure was measured via a side tap of the tracheal tube with a Validyne MP45 (±30 cmH_{2}O) pressure transducer. Tracheal flow (V˙tr) was measured with a heated screen pneumotachograph (28 mm ID) and another Validyne MP45 (±5 cmH_{2}O) transducer. All signals were lowpass filtered at a cutoff frequency of 5 Hz, then sampled at 20 Hz. Apparent lung input impedance (Zapp) at the NSND input frequencies was calculated as the ratio of the crosspower spectrum ofV˙tr and Pao and the autopower spectrum ofV˙tr.
Simulation studies.
To test the abilities of our IHNL models, we carried out the following simulation studies. The pressureflow relationship of the entire lung was simulated by using a model that includes an asymmetric branching airway system based on the Horsfield studies (15) modified to allow for inhomogeneous changes in airway diameters (20). The airway tree terminates in alveolar tissue elements that were also modified to incorporate nonlinearity. Briefly, the alveolar tissue model contains an alveolar compartment that is the parallel combination of the constantphase tissue model (Eq. 3 ) and alveolar gas compressibility. The tissue properties are uniformly distributed over all terminal airway elements. The total input impedance of the model is then determined by combining the impedance of the alveolar tissue elements and the airways in the proper parallel and serial fashion (8). The model also accounts for airway wall viscous, elastic, and inertial properties as described previously (10). Inhomogeneous peripheral airway disease was simulated by decreasing the airway diameters in various fractions of the last five generations of the peripheral airways. We chose to constrict all airways subtended by order 46 in the dog lung (20). This simulated a lung disease with 72% of the airways constricted and 28% kept at the baseline diameters. The constriction levels were set to 30, 40, 50, and 60% decrease in diameters from their baseline values at functional residual capacity (FRC). After the input impedances were obtained, the time domain pressure output of this linear model was calculated and Eq. 4 was used to add nonlinearity. This is equivalent to replacing the parallel set of airway pathways terminated by the constant phase model in Fig. 2 with the above Horsfield structure. Zapp of the nonlinear model at the NSND frequencies were then calculated as described below. These data were fit with our HNL and IHNL models of the previous section and used the modelfitting approach described below. The tissue parameters (G, H,a _{2},a _{3}) were kept constant in the simulations and were evaluated as a function of airway constriction as obtained from the fittings.
Computations.
According to Fig. 2, the airway opening flowV˙(t) is the input signal to the model. The predicted output pressure of the model, Pao,m(t), is calculated as follows. TheV˙(t) is subject to the linear subsystem to obtain a temporary pressure output, P_{Lin}(t). First, P_{Lin}(ω) =V˙(ω)Z_{L}(ω) is calculated, and then P_{Lin}(t) is obtained by taking the inverse Fourier transform of P_{Lin}(ω). The input impedance of the linear subsystem Z_{L}is defined in Eqs. 1012 . The P_{Lin}(t) signal is then passed through the nonlinear subsystem with a thirdorder polynomial in Eq. 4 to obtain the total pressure drop across the model, Pao,m(t). The apparent input impedance for both fitting measured data and simulating data (i.e., inverse and forward models, respectively) Zapp,m was also calculated at the NSND input frequencies (i.e., the frequencies at which the input flow contained energy).
Modelfitting approach.
Because of the design of the NSND frequencies, nonlinearity in Zapp at the NSND frequencies is minimized (32). Nevertheless, to identify nonlinearities from an OVWNSND input flow signal, we can fit the time domain pressure. The time domain recording does retain the effects of all nonlinear harmonic influence from which the the nonlinear coefficients can be identified (33). During lung disease, however, the dominant change in Z_{L} is the increased frequency dependence of Rl and El, which can be because changes in tissue properties and/or a consequence of increased airway constriction inhomogeneity. A twocompartment linear impedance model mimicking inhomogenity has been found to provide much better fit to lung input impedance data than the singlecompartment model during constriction in frequency domain (20). This suggests that inhomogeneity may be better identified in frequency domain. Thus we hypothesized that, to identify both inhomogeneity and nonlinearity, we must simultaneously fit our models to both time domain pressure data and apparent input impedance at the NSND frequencies. Accordingly, the optimal parameters of the model are obtained by minimizing the following mean square error
RESULTS
Validation of combined timefrequency domain modelfitting approach.
Consider, first, the fit of linear and nonlinear forms of the tissue models with and without airway inhomogeneities only to the time domain pressureflow data (i.e.,W_{t} = 1) from a dog during peak bronchoconstriction. The peak response data showed the largest degree of inhomogeneities (22). The model fitting errors of the four models are shown in Fig.3 A and clearly indicate that including nonlinearities was far more important than including inhomogeneities for improving the quality of fit (i.e, reducing the residuals vs. time, Fig.3 A). Compared with the linear model, the error was ∼50% lower when only nonlinearities are accounted for (HL and IHL vs. HNL or IHNL) and insignificantly lower when only inhomogeneities are included (HL vs. IHL and HNL vs. IHNL). This occurs because the time domain data actually deemphasizes the influence of inhomogeneities. Overall, the mean error decreased by only 1.6% when inhomogeneities were added to the nonlinear model. Furthermore, there were no significant improvements in the predicted Rl and El when inhomogeneities were included in the model (Fig. 3, B andC). In fact, the prediction of the Rl spectra are nearly identical with all models but, more importantly, no model could capture the rate of increase of El. In summary, in the modeling of time domain data, it is essential to include nonlinearity in the model; however, inhomogeneities cannot be identified from time domain data only.
We next examined whether inhomogeneities can be more reliably identified by including the frequency domain data in the identification. We repeated the above model fitting by using the combined domainfitting technique (i.e., Eq.13 ) withW_{t} = 0.1 andW_{f} = 0.9. The nonlinear coefficientsa _{2} anda _{3} obtained from the combined domain fitting (−0.0395 and −0.00138, respectively) were very close to those obtained from the time domain fitting (−0.0364 and −0.00113, respectively). Figure4 demonstrates that now the prediction of the Rl and El spectra is excellent but only with the inhomogeneous models, and this is true regardless of whether nonlinearities were included. It is worth appreciating that fitting a nonlinear model only in the frequency domain at selected NSND frequencies (26) will result in biased parameters. The reason is that, at the NSND frequencies, the higher order harmonic distortions and crosstalk are minimized and only certain crosstalk terms contribute to the data (32). Hence, the estimated nonlinear coefficients may not capture the essence of the nonlinearity. Overall, we can draw the following conclusions. When OVWNSND flow input is applied at the airway opening, identification of both airway inhomogeneities and tissue nonlinearities requires the combined timefrequency domain fitting approach.
Model form and parameters.
Among all models (the inhomogeneous models with the three different airway resistance distributions), the inhomogeneous model with a hyperbolic resistance distribution (μ = 1) had the best fit to the peak response as well as in control and the late response. The inhomogeneous models are always superior to the homogeneous models. The errors of the different models are compared in Table1. In peak response, and even in control, there are large error reductions for the inhomogeneous models over the homogeneous model (>30%). Moreover, Fratio tests show that the IHNL (with μ = 1) model is statistically superior to the HL model for all four dogs in control, peak, and late response, and it is also statistically superior to the HNL model for all dogs in peak response, three dogs in late response, and only one dog in control. Thus, while the inclusion of nonlinearity was essential for all conditions, inhomogeneity was needed mostly in the peak and late responses.
The parameters obtained from fitting the best model, i.e., IHNL withn(R) hyperbolically decreasing (μ = 1), to all dog data are given in Table2. Between peak response and control, both G and H significantly increased by 209% (P< 0.001) and 63% (P < 0.01), respectively. The Iaw slightly and statistically nonsignificantly decreased by 20%. The expected value of R, E[R], increased from control to peak response in all dogs and on average by 419%. The range between R_{a} and R_{b} should be related to the degree of airway inhomogeneities. Indeed, the widest range occurred at peak response, which was an ∼500% increase compared with the control range. The nonlinearity coefficientsa _{2} anda _{3} decreased from control to peak response by 40 and 50%, respectively. In late response, the parameters generally tended to become closer to their control values, and only H remained statistically significantly higher than its control value.
Simulation study.
Simulated pressureflow and impedance data (both control and constricted) obtained from the Horsfield airway structure terminated with the nonlinear tissue model were fit with our HNL and IHNL models. The errors obtained by the HNL and IHNL models are shown in Fig.5. In control, both model provided an excellent fit with practically no systematic error. For constrictions >30%, the reduction in error due to including inhomogeneities in the model was quite significant. The estimated parameters from both models are shown in Fig. 6. Recall that the tissue parameters were not altered in the simulations. The tissue damping G (which was set to be 3 cmH_{2}O/l in the simulations) was recovered by fitting the IHNL model within 30% at all constriction levels. Alternatively, G was substantially overestimated (>200%) by HNL model at 60% constriction. Tissue elastance H and nonlinear tissue coefficienta _{2}, which were set to 18 cmH_{2}O/l and 0.05 cmH_{2}O, respectively, were accurately estimated up to 50% constriction, and at 60% constriction the errors were still within 10–15% when the IHNL model was fit. The coefficienta _{3} (0.005 cmH_{2}O^{−2}) was overestimated by both models but to a much lesser extent by the IHNL model (errors 60 vs. 90%). The E[R] in the IHNL model increased much more than parameter R in the HNL model. The Iaw parameter decreased in both models but to a lesser extent when the IHNL model was fit to the simulated data.
DISCUSSION
The primary purpose of this study was to develop analytic tools that can better distinguish between real tissue changes and virtual tissue changes due to heterogeneous airway constriction after an agonistinduced bronchochallenge. The basis for our approach was the finding that frequency and amplitude dependences of Rti and Eti are coupled during challenge (33). Additionally, we also hypothesized that airway inhomogeneities contribute to the frequency dependence of Rl and El and that this is primarily a linear phenomenon. Thus models that simultaneously account for nonlinear tissue viscoelasticity and airway inhomogeneities could have the potential to separate the contributions of tissues and airway inhomogeneities to Rl and El. Our findings indicate that in bronchochallenged dogs, although inhomogeneities are significant in both peak and late responses, there are also real changes in lung tissue properties. Before discussing the physiological implications, we must first examine to what extent our model selection is appropriate when fitting measured data and what the limitations are of the our modeling approach on the basis of analyzing simulated data.
Model selection.
With regard to the inhomogeneous airway compartment, we note that our purpose was to develop a reasonably simple model that still captures the distributed parallel nature of inhomogeneities. We only dealt with small integer values of the exponent μ in the airway resistance distribution (Eq. 8 ) because this led to closedform solutions of the probabilistic airway compartment. The resultant model complexity thus allowed iterative minimization of the error function in Eq. 13 , while maintaining identifiability with respect to available data. Indeed, in most cases the global optimization (6) provided a single unique set of parameters. More importantly, we hypothesized that it is the alterations in the distribution of pathway resistances that contributed most to the increased inhomogeneities during constriction. This is a key point because the underlying assumption of all data analysis is that nonlinearity comes exclusively from the tissues and that airways only contribute to time constant inequalities. To test this, we included a Rohrertype airway pressureflow nonlinearity in the IHNL model. Similar to our previous study (33), this added complexity did not significantly improve the fits in any of the conditions and often resulted in nonuniqueness. This finding thereby supports our key assumption that inhomogeneities contribute to the data primarily as a linear phenomenon.
In our previous study, we found that, among six different nonlinear block structures (e.g., a Hammerstein structure, which is an inverse cascade connection of the LN structure in Fig. 1, i.e., NL), a Wiener structure for the lung tissues in series with a linear singlecompartment airway model was the most appropriate for describing the pressureflow behavior in tidalventilated dog lungs in control and late response to histamine infusion (33). In this study, we combined the previous model with a distributed airway compartment. This IHNL model provided a significantly improved fit to the data in both the peak and late responses compared with the HNL model. To verify that the best lung model was indeed a Wienertype structure, we also fit the Hammerstein structure with and without airway inhomogeneities. Although not shown, these additional modeling efforts during both peak and late responses remained inferior to the Wiener structure regardless of the airway compartment included in the model (i.e., singlecompartment or distributed model). Indeed, on average, the errors by using the Wiener structure were 45% smaller than the errors with the Hammerstein structure.
The linear part of the lung tissue model was the constantphase model given by Eq. 3 . This model has been shown in several previous studies to provide excellent fits to small and tidallike amplitude tissue impedance data (3,11–13, 20–22, 27, 31, 33) as well as pressure and stress relaxation data (3, 31). To verify that this was the case in the peak response, we replaced the constantphase model by a Kelvin model, which did not improve the fits despite the increased number of parameters (i.e., 2 vs. 3). Furthermore, on the basis of our previous study (33), here we combined the various airway compartments with only a thirdorder zeromemory nonlinear block (N block in Fig. 1). As an additional check on the model validity, we also verified that the thirdorder polynomial nonlinearity was still sufficient. Although an arctangent form or a secondorder polynomial nonlinearity gave worse fits, including a fourthorder polynomial did not improve the fits statistically significantly.
Limitations of the modeling approach.
We established the applicability of the modelfitting approach by examining its ability to recover the tissue parameters (G, H,a _{2} anda _{3}) from simulated data in which the tissue parameters were kept constant and inhomogeneous airway constriction was gradually increased. The model used in the simulations is an anatomically consistent structural model based on Horsfield’s morphological data (15). In several previous studies, we extended this model to be able to simulate inhomogeneous airway constriction (20). By adding nonlinear features to lung tissue properties, we have created, to our knowledge, the most faithful description of oscillatory pressureflow relationships in the lung. Our best model with a hyperbolic distribution of parallel airway resistances is able to recover all tissue properties within 15% error when the constriction in 72% of the peripheral airway diameters in the last five generations is <60%. At 60% constriction, most parameters are still reasonably recovered, the largest errors occurring in G anda _{3}, which are overestimated by 11 and 57%, respectively. Note, however, that the HNL model provided much higher errors in G anda _{3} (194 and 104%, respectively).
At constriction levels above a 60% diameter reduction, the situation becomes worse, and other parameters (e.g., H and a _{2}) will also become biased. The reason for the increased error in the estimated parameters is due to the fact that as peripheral airway constriction increases the impedance of the peripheral airways increases significantly. When peripheral airway impedance becomes comparable to central airway wall impedance, less flow enters the periphery and more flow is shunted into airway wall motion. Indeed, such a mechanism has been proposed to account for lung impedance data during severe bronchoconstriction (2, 24).
Another limitation of the model analysis is the sharp serial distinction between airways and tissues. It is conceivable that when small airways constrict they will distort the surrounding parenchyma and that, as a consequence, tissue mechanical properties will change. Although correlations between airway resistance and tissue resistance or compliance during challenge have been reported (19), to our knowledge, such an airwayparenchymal interaction has not been implemented in a model of the oscillatory pressureflow relationship of the lung. We examined this possibility by calculating correlations between the changes in the airway and tissue parameters in our model. No significant correlations were found between airway and tissue parameters. However, G and H showed some correlation (correlation coefficient = 0.63). This indicates that because this model takes into account the heterogeneities of pathway resistances, the tissue parameters are correlated with each other more than with the airway parameters. Correlation between G and H conforms with the structural damping hypothesis that postulates that dissipative and elastic phenomena are related at the elementary level (9).
Physiological implications.
Our IHNL model produces much smaller values for tissue damping G and and slightly smaller values for tissue elastance H than the HL model fit to the same data in the frequency domain (22). For example, in peak response, the IHNL model produced reduced estimates of G and H compared with the HL model by 44 and 20%, respectively. We attribute this to the fact that the inhomogeneous model can account for artificial increases in tissue resistance and elastance. Furthermore, for both control and late responses, the values of G obtained from the IHNL model were smaller that those obtained from the HNL model in our previous modeling study (33): the percent decrease in G was 14.6% in control and 25.6% in late response. The greater decrease in G is consistent with a stronger influence of inhomogeneities in late response. Additionally, both the HNL and the IHNL models in this study produced very close values for H as well as fora _{2} anda _{3}, with their values in good agreement with those found in our earlier study in which the data were only fit in the time domain (33). The implications are that the hysteresivity coefficient, η = G/H, is also reduced by ∼20% compared with that predicted by the HL or HNL models in control, late and peak responses.
By using the IHNL model, G and H in the late response were elevated compared with their control values (Table 2). Although the magnitudes of these elevations were similar to those previously found (33), these increases reached a statistically significant level for only H. This indicates that in our previous study (33), part of the increase in G in the late response may have been attributed to inhomogeneities. Additionally, parameters a _{2} anda _{3} slightly decreased in the late response but again not statistically significantly. Because these results were obtained by a model that accounts for inhomogeneities, they suggest the existence of real tissue changes in the late response to histamine infusion. Thus, consistent with our previous finding, we conclude that increased nonlinearity in late response after histamine challenge as measured by the Vt dependence of apparent impedance is not due to changes in tissue nonlinearities but to an increase in the magnitude of the linear tissue impedance (33). This is also supported by the fact that H did not correlate with the nonlinearity parametersa _{2} anda _{3}.
During peak response, both G and H increased statistically significantly compared with their control values (Table 2). The nonlinear coefficients decreased again, but their decrease did not reach a statistically significant level. It is possible that part of the increases in G and H was still a result of inhomogeneous airway constriction or central airway wall shunting. However, G obtained from the IHNL model decreased considerably by 26% compared with that of the HNL model. In addition, the IHNL model was able to recover G and H within 15% when the Horsfield structure undergoes a 60% inhomogeneous peripheral constriction (Fig. 6,A andB). Thus this presents indirect evidence that in dogs real tissue changes do occur at the level of the parenchyma in both peak and late responses during bronchochallenge. This is in agreement with some recent findings (13, 23). Nevertheless, our study suggests that the tissue response may be considerably less than previously thought (1,18, 19, 22, 26, 28).
The airway structure in our model consists of a set of parallel pathways. The most appropriate distribution function for the parallel airway resistances was a hyperbolically decreasing distribution with a long tail toward the higher resistances. Interestingly, this implies that airway conductance is uniformly distributed. A uniform distribution of airway conductances has been used to predict the effects of peripheral inhomogeneities on input impedance by Hantos et al. (12). The expected value of the distribution E can be quite different from the arithmetic mean of its left and right limits, i.e., (R_{a} + R_{b})/2. As constriction level increases in the simulated data (Fig. 6), E becomes increasingly higher than the parameter R from the HNL model. One might expect that the width of the distribution of resistances in this model correlates with the degree of inhomogeneities. The simulation results show that by using the Horsfield structure for the control case, the estimated width of the distribution is negligible and hence the model predicts a fairly homogeneous system (Fig.7 A). When constriction is increased in 72% of the peripheral airways of the Horsfield structure, the width of the estimated airway resistance distribution also increases. This extended tail of the distribution increases the frequency dependence of both Rl and El at low and higher frequencies, respectively (see Fig. 4). Hence the IHNL model can partition Rl differently from any homogeneous model, and the error in G is much smaller than that provided by the HNL model (Fig.6 A). A similar tendency can be seen from the fits to the measured data: both the expected value (Table 2) and the width of the estimated distribution functions increase from control to peak response and then decrease from peak to late response (Fig. 7 B). Experimentally, the quantification of airway inhomogeneities seems quite arbitrary by using alveolar capsules because of the inevitable undersampling of alveolar pressures. On the other hand, our approach provides an average behavior of the airwaytissue structure, and thus the resistance distribution should reflect changes in peripheral time constant inequalities. Indeed, the IHNL model seems to be sensitive to inhomogeneous peripheral constriction simulated by the Horsfield structure.
In summary, we presented an approach to distinguish between lung tissue changes and airway inhomogeneities by using a model that features nonlinear tissue viscoelasticity as well as airway inhomogeneities. For this approach to provide sensible results, we had to develop a combined time and frequency domainfitting procedure because time domain data are dominated by nonlinearity, whereas fitting in frequency domain is necessary for the estimation of the contribution of airway inhomogeneities. Another method to distinguish tissue changes from airway inhomogeneities would be to use gases of different viscosities before and after steadystate bronchoconstriction (21). However, this presents a substantial increase in experimental and protocol complexity and may certainly prove difficult in asthmatic subjects. On the other hand, our approach is based on a single measurement made at the airway opening and does not require the use of the invasive alveolar capsule technique or different gases. Therefore, this technique offers great potential for use in human subjects in dynamic challenge studies. Nevertheless, future studies should examine and possibly include additional mechanisms such as tissue inhomogeneities, airway tissue interdependence, and central airway wall shunting.
Acknowledgments
This work was supported by National Science Foundation Grants BCS9309426 and BES9503008 and National Heart, Lung, and Blood Institute Grant HL50515.
Footnotes

Address for reprint requests: B. Suki, Dept. of Biomedical Engineering, Boston Univ., 44 Cummington St., Boston, MA 02215 (Email:bs{at}enga.bu.edu).
 Copyright © 1997 the American Physiological Society