Abstract
The maximal expiratory flowvolume (MEFV) maneuver is a commonly used test of lung function. More detailed interpretation than is currently available might be useful to understand disease better. We propose that a previously published computational model (Lambert RK, Wilson TA, Hyatt RE, and Rodarte JR. J Appl Physiol 52: 44–56, 1982) can be used to deduce, from the MEFV curve, the serial distribution of airway areas in the larger airways. An automated procedure based on the simulated annealing technique was developed. It was tested with modelgenerated flow data in which airway areas were reduced one generation at a time. The procedure accurately located the constriction and predicted its size within narrow bounds when the constriction was in the six most central generations of airways. More peripheral constrictions were detected but were not precisely located, nor were their sizes accurately evaluated. Airway areas of generations upstream of the constriction were usually overestimated. The procedure was applied to spirometric data obtained from eight volunteers (4 asthmatic and 4 normal subjects) at baseline and after methacholine challenge. The predicted areas show individual differences both in absolute values, and in relative distribution of areas. This result shows that detailed information can be obtained from the MEFV curve through the use of a model. However, this initial model, which lacks airway smooth muscle, needs further refinement.
 flowvolume
 human
 spirometry
 asthma
 maximal expiratory flowvolume curve
since its introduction more than 40 years ago, the maximal expiratory flowvolume (MEFV) curve has become the most widely used test of lung function (9, 11). However, despite our present understanding of the mechanisms underlying the maximal flow phenomenon (6, 7, 9, 10, 12, 14, 15, 17, 19, 20, 25, 27), interpretation of the MEFV test is often vague and uncertain. We hypothesized that it may be possible to obtain more precise information about airway properties from the test. As with most test interpretation, this must be done in terms of a model. We chose to use an established computational model for forced expiratory flow (19).
This model for forced expiration from humans (19) has been used previously to examine the effects of various airway abnormalities on the MEFV curve (14, 15). The modeling was implemented by manipulating the input parameters to the model to produce changes in the MEFV curve. This procedure yielded insights into the link between the mechanical and geometric properties of airways and lung function.
An attempt at solving the inverse problem (What airway properties would give rise to a given MEFV curve?) was unsuccessful (16). There were two reasons for this. The first was the application of a linear approach to a highly nonlinear problem. The second was that the number of points used in the MEFV modeling were far too few to find optimal values for the model parameters. A much more recent study of the inverse problem has shown that a model for forced expiration from infants can be applied by manually manipulating parameters to predict individual flowvolume curves (18). We reasoned that with the huge increase in computing power in recent years it should be possible to automate the procedure.
In this paper, we present an automated approach to the inverse problem that yields information on the serial distribution of airway lumen area. We start by presenting data that indicate how modelpredicted maximal flow is affected by large excursions in the values of the input parameters. We then apply this information to the design of a software routine that adjusts maximal airway area generation by generation until an optimal match is found between model prediction and the data. We illustrate its use both with modelgenerated data and with data obtained from human subjects.
REVIEW OF THE COMPUTATIONAL MODEL
The model we used is that developed more than 20 years ago by Lambert and colleagues (19). It contains 17 generations of dichotomously branching airways based on Weibel's model A (31). The mechanical behavior of each generation is specified by normalized areapressure (αP) curves, with α defined as A_{i}/A_{i}. A_{i} is the lumen crosssectional area, and its maximal value is indicated by the asterisk (3). A value of A_{i} is assigned to each generation. Each curve is described by two hyperbolas matched for intercept and slope at P = 0 where P is the transmural pressure difference. α_{0} is the value of α at P = 0. P_{1} and P_{2} are the pressure asymptotes of the two hyperbolas and are inversely proportional to the slope of the curves at P = 0 (α_{0}), and n_{1} and n_{2} are shapechanging constants. Additionally, the length is specified for each generation. A similar pair of equations together with a value for vital capacity (VC) is used to specify the model's pressurevolume curve. In all, 85 parameters are used to specify the areapressure characteristics of the 17 generations of airways (Table 1).
The fluid mechanics of this model are based on an expression for the pressure gradient (dP/dx) along a compliant airway (Eq. 1). (1)
In Eq. 1, f indicates the dissipative pressure loss per centimeter of tube length, and U and c are the local values of flow speed and wave speed, respectively. The expression used for f is that developed by Reynolds (26) for a human lung (Eq. 2). (2) V̇ is the volume flow in an airway, μ is the viscosity of the gas, and Re is the Reynolds number of the flow in that airway. Eq. 1 is integrated along an airway to calculate the total pressure drop along that airway generation. Because all airways in a generation are identical and in parallel, this is the same as the pressure drop across that entire generation of airways. To calculate the pressure change caused by merging flows in a junction, the Bernoulli equation is solved for the pressure in the downstream airway. It is assumed that the pressure dissipation in the junction is accounted for by f in Eq. 2. This provides a starting pressure for the integration along the next generation of airways. The procedure is repeated along each Weibel generation of airways (31) in turn, starting from alveolar pressure at the upstream end of generation 16, the most peripheral nonalveolated generation, until the downstream end of the intrathoracic trachea is reached. Computation is stopped by either U rising to within 0.2% of the local value of c or the pressure in the airway falling >100 cmH_{2}O below pleural pressure. Under these conditions, and in all cases in which n_{1} > 0.4 in the flowlimiting generation, each isovolume pressureflow curve is on its plateau (15).
SENSITIVITY OF THE MODEL TO ITS PARAMETERS
The sensitivity of the model to the most important of its parameters was examined by altering each of the parameters in turn, generation by generation, and calculating the new predicted maximal flow in each case. The ratio of the new predicted flow to the baseline flow was calculated and taken as the output variable for this part of the study. The four parameters studied were A^{*}_{i}, α_{0}, α^{′}_{0}, and bronchial length (BL). The earlier study clearly showed a lack of sensitivity to n_{1} and n_{2} (16).
The baseline value of each parameter was multiplied by factors of 0.33, 0.50, 0.67, 1.5, 2.0, and 3.0. The results showed a gradual evolution in sensitivity with change in parameter. As in the previous study, it became apparent that the model is markedly less sensitive to the values of α^{′}_{0} and BL than to A^{*}_{i} and α_{0}. Therefore, only the results for A^{*}_{i} and α_{0} are shown here for factors of 0.33 and 3.0 (Figs. 1 and 2). Multiplying by a factor of 0.33 in any airway generation causes large changes in maximal flow but over only a small range of lung pressure (Pl) (Figs. 1, A and B, and 2, A and B). On the other hand, the effect of increasing either of these parameters by a factor of 3.0 makes only a small difference to the predicted maximal flow (Figs. 1, C and D, and 2, C and D). For instance, increasing the value of A^{*}_{i} by a factor of 3 has a significant effect only in the most peripheral airways (Fig. 1) at low lung inflations. Even then the effect is small, with the greatest change being an increase in maximal flow to 20% above baseline when the A^{*}_{i} of generation 12 is increased by a factor of 3 at a Pl of 0.13 cmH_{2}O. When the same parameter is decreased by a factor of 3 at the same Pl, flow is reduced by 65% of baseline.
COMPUTATIONAL SCHEME
We used the technique known as simulated annealing as the basis of a software package (FVfits) with which to search for the optimal fit between the model predictions and a real MEFV curve by adjusting model input parameters.
Our procedure starts by using the model to generate an MEFV curve from an initial guess at the required parameter set, and from this curve a cost function that is a measure of the goodnessoffit of this model result to the data is computed. A new parameter set is then generated from the old by multiplying each of the parameters by its own randomly generated factor. A new MEFV curve and cost function are then calculated. Whether the perturbation to the parameter set is acceptable is determined probabilistically in analogy with the thermal processes involved in annealing of metals. In annealing, a hot metal is slowly cooled to achieve the lowest energy state. The probability of a transition from a state with energy E_{1} to a state with energy E_{2} is governed by the Boltzmann factor, exp[−(E_{2} − E_{1})/kT], in which T is the absolute temperature and k is the Boltzmann constant. As the system cools, the probability of the system making a transition to a highenergy state becomes less. Slow cooling allows thermal processes to explore a vast number of configurations to find that with the minimum energy. The analogous computational procedure (13, 22) was adapted for use in this study (24). It uses the cost function as the measure of system “energy” that is to be minimized as well as the basis for an initial value of the analog of temperature that is to be reduced as the process proceeds. Increases and decreases in cost function are accepted in a probabilistic manner. Large increases in cost function are less likely to be accepted than smaller increases. In this way, the parameter space is explored to find (one hopes!) the global minimum of the cost function. The scheme for accepting or rejecting changes and for lowering the temperature as the minimum value of the cost function is approached is known as the annealing schedule. The schedule is unique to each problem, and the best schedule has to be found by trial and error.
Our cost function is formed by calculating the sum of the absolute differences between predicted flows and data. We found that assigning an initial temperature of 2% of the initial cost function worked satisfactorily for our problem. Following Metropolis and colleagues (22), we accept all decreases in cost function and apply the probabilistic acceptance only to increases. We decrease the “temperature” in steps of 10% of the current temperature (that is, the new temperature is 90% of the previous temperature) until either we have tried 200 parameter sets without finding a further decrease in cost function or the temperature falls to less than a preset limit.
The data we planned to use had flow at 2% intervals of forced vital capacity (FVC). In adults, expiration is flow limited from ∼15% VC expired until the end of expiration. We use the final 45 data pairs (volume and flow) from an expiration as our data set and use the 17 values of A^{*}_{i} as our parameter set. All other model parameters (Table 1) are held constant at baseline values. (However, BL is different for different lungs because it is scaled on the cube root of lung volume.) Thus we have 2.65 data points/parameter value. This is a light constraint on the values of A^{*}_{i} (think of fitting a straight line to 5 points).
Annealing starts with an initial guess at the A^{*}_{i} distribution in which each value is greater than the expected optimized value in that generation. During optimization, for each new value of A^{*}_{i} to be accepted, it must lie between the initial value and 10% of that value for that generation.
Our MEFV model has no time dependence in it. Thus the maximal flow at any lung volume can be calculated without knowledge of the flows at greater volumes. This is useful to us in this study for the following reason. The peripheral airways are the most important determinants of maximal flow at small lung inflations (Figs. 1 and 2). These airways also influence flows at all volumes. Thus if the values of A^{*}_{i} in the peripheral airways are not right, then no other values can be right. These areas should, therefore, be estimated first by examining only the last part of an expiration. Hence the model is initially run between Pl = 0 and 1.5 cmH_{2}O with only the values for A^{*}_{i} in generations 8–16 being manipulated. Continuing this approach, optimal values of areas for generations 4–7 are then found by using data in the Pl range 1.5 to 6 cmH_{2}O and leaving the values for generations 8–16 at their optimal values. The model is finally run over the total 87% VC expired range to find optimal values for generations 0–5 while maintaining the existing optimized values for the other generations. Including generations 4 and 5 in two parts of the procedure is done because the results shown in Fig. 1 indicate that these generations exert their maximal influence from ∼3 to 9 cmH_{2}O.
The optimized fit depends on the value used to start the random number procedure (the seed). Therefore, for each MEFV curve used, the process is repeated 20 times using 20 different seeds to provide a measure of the variability of each fitted value of A^{*}_{i} (4).
APPLICATION OF FVfits TO SYNTHETIC DATA
We tested FVfits by using synthetic data generated from the same computational model used in the fitting procedure. Thus, at least in theory, it should have been possible to obtain a perfect fit to these data. We generated eight test curves by reducing the baseline value of A^{*}_{i} by 50% in each of generations 0–7 and then ran the model to obtain a synthetic MEFV curve with each (known) airway abnormality. We ran the model at the same 2%VC increments that were present in the real data. These synthetic flow data were then used as input to FVfits. In each case, 20 attempts were made to find an optimal fit to the input MEFV data by manipulating A^{*}_{i} values from an initial set in which each value of A^{*}_{i} was twice that in the baseline model. The MEFV curves found by using FVfits were excellent matches to the input curves. The mean of the absolute residuals for the best fit (defined as that fit of the set of 20 that had the smallest cost function) curves ranged from 3 ml/s (constricted generation 5) to 12 ml/s (constricted generation 7). The constricted generation was clearly identified in generations 0–5. The results for generations 6 and 7 were more ambiguous.
To make the test of the software more realistic, a small amount of normally distributed noise (mean 0 ml/s, SD 30 ml/s) was added to the synthetic flow data, and the optimization procedure was repeated as described above. The addition of noise resulted in greater parameter variability. However, the constricted generation was still clearly identified in generations 0–5. Parameter variability for the simulations with added noise is summarized in Fig. 3. Results for generations 0, 2, 3, and 4 were similar to those shown for generations 1 and 5. It is clear that not only was the constriction identified, the area at the constriction was confined accurately within a narrow range for generations 0–5.
The MEFV curve contains no information about airways downstream of the choke point. This point moves toward the periphery as the lung deflates. Thus the most central location of the choke point at the greatest lung inflation defines the most central airway for which information can be deduced. Because of this, Fig. 3 shows no A^{*}_{i} predictions for airways downstream of the most central choke point. For the constriction in generation 7, this point was downstream of generation 7 and thus there are A^{*}_{i} predictions for two airway generations downstream of generation 7. For the constriction in generation 1 there is no tracheal information. When the constriction was in generation 5, some fits had the most central choke point downstream of generation 5, but most did not. Thus the plotted predictions for A^{*}_{i} do not contain results for generations more central than generation 5.
We constricted (by 50%) both generations 5 and 6 together as well as generations 7 and 8 together and added noise to the flow data. In the first case, the fitting procedure put almost all of the constriction in generation 5, whereas in the second case the constriction was spread over generations 6–9 (as indicated by the median area values) rather than being restricted to generations 7 and 8 (Fig. 3). Constrictions peripheral to generation 8 were not localized accurately by our procedure, but the presence of a constriction was apparent.
Whereas FVfits appears to find the constriction with good accuracy in generations 0–5, it does not provide accurate assessments of areas upstream of the constriction. In almost all cases the median values of A^{*}_{i} are greater than the input values (Fig. 3).
EXPERIMENTAL PROTOCOL FOR HUMAN SUBJECTS
The experimental study was performed at Mayo Clinic, Rochester, MN under approval of the Mayo Institutional Review Board. Data from the normal subjects from this study have been presented (1). Four normal subjects and four subjects with a history of asthma participated after signing the approved informed consent form. Subjects were asked to swallow an esophageal balloon, which was used in measuring pulmonary resistance (Rl) at various times during the study by the Mead and Whittenberger method (21) between flow limits of ±2 l/s, as described by our laboratory (5). A methacholine (MCh) challenge study was performed on all subjects as follows. First, prechallenge spirometry was performed using custom software and a screentype pneumotachograph precalibrated by using a 3liter gasfilled syringe, as described (5). At least three maneuvers were performed, and the maneuver with the best FVC and forced expiratory volume in 1 s (FEV_{1}) was chosen to represent baseline data.
Subjects were then asked to breathe mist from a nebulizer for 1 min, while following the beat of a metronome set to 10 breaths/min. Three minutes after stopping the mist, we recorded normal breathing, including esophageal balloon measurements to calculate Rl postMCh, followed by an FVC maneuver at the 4th min. After subjects relaxed for 3 more min, we again recorded normal breathing with esophageal manometry followed by a final FVC maneuver. For this report, we analyzed the MEFV curves from both the prechallenge and the first postchallenge maneuvers.
This procedure was repeated with the nebulizer primed first with 3 ml of normal saline, and then with increasing doses of MCh: 0.1, 0.4, 1.6, 6.4, and 25 mg/ml. Challenge testing was stopped when the subject's FEV_{1} in the 8th min postMCh administration fell below 80% of the preMCh test. Three of the four asthmatic subjects experienced >20% drop in FEV_{1} at a dose lower than 25 mg/ml, and one of the normal subjects experienced a >20% drop at the 25 mg/ml dose.
EXPERIMENTAL RESULTS
PreMCh data for all subjects are shown in Table 2. The spirometry and Rl data for the highest dose of MCh administered to each subject are shown in Table 3. Note normal to mildly abnormal spirometry preMCh in all of the asthmatic subjects, and completely normal spirometry in all of the normal subjects. After MCh, Rl increased dramatically in both the normal and asthmatic subjects.
In the application of FVfits to experimental data, results from two normal subjects and one asthmatic subject are shown in Fig. 4. Figure 4A shows the data and bestfit MEFV curves, and Fig. 4B shows the median area distributions underlying the model flow results. Differences in VC were accommodated by dimensionalizing the model's normalized MEFV curve by using the FVC recorded during the maneuver of interest. Results showing the effects of MCh challenge in one normal subject are shown in Fig. 5. Note in Fig. 5B that a narrowing in medium generation airways seemed to account for the nearly 50% drop in flows shown in Fig. 5A. Mean data for area distributions after MCh for normal and asthmatic subjects are shown in Fig. 6. In all cases, the model placed the constriction in generation 3 and higher for the asthmatic subjects and in generation 4 and higher for the normal subjects.
DISCUSSION
We have shown that the computational model for forced expiration can be fitted to individual MEFV curves by choosing an appropriate value for maximal lumen area in each generation of airways. We tested our procedure by generating model curves to be fitted in which the airway tree contained various known constrictions. We added noise to the flow data. When the constriction was in the central and middle airways, our procedure recovered the location of the constriction and the area at the constriction with considerable accuracy (Fig. 3). Area information for more peripheral constrictions was less accurate, as shown by the generation 7 results (Fig. 3).
Neither the areas nor the locations of twogeneration constrictions in the middle airways were accurately recovered from the synthetic flow data (Fig. 3). However, the constricted generation 5 and 6 results (not shown) gave a more central predicted constriction, predominantly in generation 5, than that given by constriction of generations 7 and 8 (Fig. 3). That is, the two constrictions were correctly located with respect to each other.
When the curve to be fitted was obtained from routine spirometry, our procedure produced a good fit to the flow data. The fitted airway area distributions show individual differences that are consistent with our understanding of the MEFV curve (Fig. 4). The curve from the asthmatic subject has greater flows than either of those from the two normal subjects at all volumes. This appears to result from greater areas in all generations. The two curves from the normal subjects differ from each other at higher lung volumes and this appears to be caused by normal subject 4 having larger cartilaginous airways than normal subject 3. This is consistent with the sensitivity results (Fig. 1). An apparently interesting feature of the results shown in Fig. 4 is the peak in the relative area distribution at generation 7 in all three subjects. However, this is an artifact of the choice of results in Fig. 4. Fits to the flowvolume curves of other subjects do not show this peak.
It is apparent from the area distributions in Fig. 3 that the median values are mostly overestimates of the true values. This was to be expected. The sensitivity study showed that increases in A^{*}_{i} have much less effect on maximal flow than decreases (Fig. 1). One of the reasons is that increasing the area of a generation frequently changes the location of the choke point. If the choke point is normally in the generation being manipulated, increasing the area results in the choke point moving out of that generation with little change in maximal flow. However, a reduction in area in the flowlimiting segment results in a substantial and areadependent reduction of maximal flow because the choke point stays in that generation. When the airways of concern are upstream of the choke point, a doubling of area produces a much smaller fractional increase in flow than the fractional decrease produced by a halving of the area. This is a result of the nonlinear dependence of airway resistance on airway area. Thus the increase in predicted flows that would be caused by several generations of airways being larger than the input values can be compensated by a single generation downstream of the enlarged generations being a little smaller than the input value.
Each optimization produced a different set of A^{*}_{i} values. Thus we had 20 different results, any one of which gave an MEFV curve that was a close approximation to the data. We took medians as a simple form of data reduction. Because the A^{*}_{i} values are not normally distributed, the usual summary statistics used to indicate the spread of data are not appropriate. Therefore we have presented modified box plots (Fig. 3) to give an indication of this spread.
The model has a total of 85 parameters to specify the tube law in the 17 airway generations. We had only 45 data pairs with which to perform our optimization. Thus we had to choose which parameter(s) should be used. We chose to use A^{*}_{i} because the model is sensitive to it. The model is almost equally (but differently) sensitive to α_{0} (Figs. 1 and 2), and we could have chosen to optimize this parameter but preferred to work with A^{*}_{i}. We did not optimize both parameters because that would have given us almost no optimizing power with only 1.3 data pairs/parameter; 2.6 pairs/parameter was little enough.
We partitioned our annealing schedule into three zones based on our sensitivity study for A^{*}_{i} (Fig. 1). By doing this instead of performing a global optimization, we obtained a faster procedure and prevented the algorithm from trying to trade off large fractional changes in the small airways against smaller fractional changes in the large airways. That is, our schedule prevented a good fit at high volumes coming at the expense of a poor fit at low volumes. We could, instead, have used a normalized cost function in which all parts of the MEFV curve were given equal weight and not partitioned the airway tree. However, our trials with such a function showed that it gave poorer convergence and required much longer execution times. We also tested the downhill simplex method described in Numerical Recipes (24). Again this was slower than our simple Metropolistype technique and yielded poorer convergence.
The final value of the cost function was different for each of the 20 optimizations performed on each data set. This shows that we most probably did not find a global minimum and therefore did not find the truly “best” fit to the data. It is not clear how much time would be needed to find this best fit because simulated annealing is somewhat akin to hopping down a gently sloping rockstrewn valley on a pogo stick looking for the deepest pothole. Finding the valley floor, or somewhere close to it, is not difficult in our problem, but finding the deepest pothole could take a very long time. It is probable that we could have obtained smaller cost functions by letting the routine run longer at the lower temperatures. However, our aim was to investigate possibility rather than to produce finely tuned software.
Fitting the model to the data is a computationally intensive procedure. A typical single fit involves ∼2,500 model runs. (We also tried using ∼4 times as many iterations with no better results.) Therefore, we investigated the possibility of using only five fits to reduce the computational time. There were no important differences between the two estimates of the median values of A^{*}_{i}. This was true in all cases investigated. Thus performing 20 fits serves to reduce parameter variability but yields essentially the same area estimates as 5 fits.
Our model cannot address the effects of parallel heterogeneity of airways. However, theoretical studies of forced expiration from lungs with parallel heterogeneity have failed to show any marker of such heterogeneity in the MEFV curve (17, 23, 30). Thus interpreting spirometric results with our homogeneous model may be all that can be reasonably achieved.
A potential use for software such as FVfits is to track changes in airway caliber during interventions or over time. We have intervention data and it is a simple matter to find sets of airway areas for the various interventions. However, because the intervention was to narrow airways by stimulating the development of airway smooth muscle (ASM) tone and the model that we use has no ASM, finding airway areas is somewhat fraught. Nevertheless, as a purely speculative exercise and to indicate what might be possible with a more appropriate model for forced expiration, we applied FVfits to our intervention data. Figure 5 shows the results obtained when normal subject 4 inspired 25 mg/ml MCh. The results appear to indicate that the effect of MCh in this subject was to produce a substantial reduction in the caliber of the 2 to 4mm airways.
We continued the speculative exercise by asking whether we could detect any overall differences in degree of airway narrowing between asthmatic and normal subjects. The median A^{*}_{i} value for each generation in each individual was calculated at the maximal concentration of MCh administered to that individual. These values were then normalized on the baseline median values and averaged over individuals by generation. The results are shown in Fig. 6, from which it could be tempting to deduce that there was a greater area response to MCh in the asthmatic subjects than in the normal subjects (as would be expected). This is consistent with the experimental data in that Rl values were greater and maximal flows were reduced more in the asthmatic subjects than in the normal subjects. However, such a deduction must be tempered with the knowledge of the size of the parameter variability (Fig. 3). Also, the lack of ASM in the model means that neither the changes in airway compliance produced by ASM activation nor the historydependent ASM tone generation have been properly accounted for (2, 8, 28, 29). On the other hand, the model is not nearly as sensitive to changes in compliance as it is to changes in lumen area. Thus the narrowing effect of ASM activation may be qualitatively reflected in the results in Fig. 6.
In summary, the most important outcome of this study is that we have shown that it is possible to obtain information about the serial airway area distribution in the cartilaginous airways from the MEFV curve under baseline conditions.
GRANTS
This work was supported in part by the Asthma and Respiratory Foundation of New Zealand, New Zealand Lottery Health Research Committee, National Institutes of Health Grants HL52230 and M01 RR00585, and the Mayo Foundation.
Acknowledgments
The authors are grateful for the programming assistance of Rock Cahan and conversations with Ted Wilson.
Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
 Copyright © 2004 the American Physiological Society