The effect of disease and respiration on airway shape in patients with moderate persistent asthma

Computational models of gas transport and aerosol deposition frequently utilize idealized models of bronchial tree structure, where airways are considered a network of bifurcating cylinders. However, changes in the shape of the lung during respiration affect the geometry of the airways, especially in disease conditions. In this study, the internal airway geometry was examined, concentrating on comparisons between mean lung volume (MLV) and total lung capacity (TLC). A set of High Resolution CT images were acquired during breath hold on a group of moderate persistent asthmatics at MLV and TLC after challenge with a broncho-constrictor (methacholine) and the airway trees were segmented and measured. The airway hydraulic diameter (Dh) was calculated through the use of average lumen area (Ai) and average internal perimeter (Pi) at both lung volumes and was found to be systematically higher at TLC by 13.5±9% on average, with the lower lobes displaying higher percent change in comparison to the lower lobes. The average internal diameter (Din) was evaluated to be 12.4±6.8% (MLV) and 10.8±6.3% (TLC) lower than the Dh, for all the examined bronchi, a result displaying statistical significance. Finally, the airway distensibility per bronchial segment and per generation was calculated to have an average value of 0.45±0.28, exhibiting high variability both between and within lung regions and generations. Mixed constriction/dilation patterns were recorded between the lung volumes, where a number of airways either failed to dilate or even constricted when observed at TLC. We conclude that the Dh is higher than Din, a fact that may have considerable effects on bronchial resistance or airway loss at proximal regions. Differences in caliber changes between lung regions are indicative of asthma-expression variability in the lung. However, airway distensibility at generation 3 seems to predict distensibility more distally.

In most cases the airway internal diameters (Din) [33][34][35][36], lumen areas (Ai) [33][34][35][37][38][39] or internal perimeters (Pi) [38,40] are evaluated. For computational models of gas-transport and aerosol deposition, idealized bronchial tree structures are developed and the above measurements used to define a network of bifurcating straight tubes [2,11,41,42]. However, it has been documented that this simplified models may not be accurate, particularly in the most proximal regions of the lung where the bronchi display transversal curvature [43,44]. Additionally, the shape of the cross-section of the airway is neither circular nor smooth [7,45] (Fig  1), a situation exacerbated in disease conditions [45][46][47]. For example, permanent airway remodeling in asthma causes the non-uniform swelling of the airway which, combined with increased mucous secretion, can form a complex lumen shape [45]. These effects need to be taken into account during functional modeling as they will clearly affect the flow.
A quantity that can be used to infer the potential functional influence of the internal shape of the airways is the Hydraulic Diameter (Dh), which is commonly used in fluid dynamics to characterize flow in non-circular channels [48]. Only a very few studies of this geometric feature have been performed in the past, with some concentrating in the extra-thoracic regions [49,50]. Recently, Choi et al. [38] have calculated regional Dh values and found that, when normalized to a predicted trachea diameter, Dh decreases in asthmatic bronchi when compared to healthy bronchi and correlates with certain Pulmonary Function Tests (PFTs).
In the current study, the influence of lung volume on the Dh, Ai and lumen internal perimeter Pi was examined. A set of bronchial measurements were performed in a mild intermittent or persistent group of asthma patients during methacholine (Mch)-induced broncho-constriction who were scanned with High Resolution Computed Tomography (HRCT) at two levels of lung inflation: mean lung volume (MLV) and total lung capacity (TLC). Subsequently, the Dh of individual airways was calculated using the measured average Ai and Pi and its correlation to the average Din was investigated. Finally, the distensibility of the airways was calculated and the correlation of lung volume change (ΔV/V) to lumen area change (ΔAi/Ai) was explored.

Medical image acquisition
The imaging protocol was completed at the Massachusetts General Hospital with IRB approval (Application No. 2007P000493) and is described in detail in [24]. In brief, the study included 24 volunteers diagnosed with mild intermittent or mild persistent asthma, as defined by the NIH Global Initiative for Asthma [51]. Written informed consent was provided by all volunteers. The first patient was recruited in December 2010. Demographic parameters are provided in Table 1. All subjects displayed reversible obstruction with inhaled albuterol. During an initial screening visit, the provocative concentration of MCh required to cause a 20% drop in the subject's FEV1 (PC20) was estimated. On a second visit that same concentration was administered over 5 deep breaths by a DeVilbiss nebulizer and Rosenthal dosimeter (model 646, DeVilbiss Healthcare, Somerset, PA) with the subject already placed in a supine position within a PET/CT scanner (Biograph 64, Siemens AG). Two HRCT images were then obtained during breath hold; one at mean lung volume (MLV) and one at total lung capacity (TLC), respectively 5 and 30 minutes after the MCh challenge. The scanner was used in a helical mode to acquire 64 slices per rotation with a 0.6mm  collimation. The imaging parameters were 120kvp, 80mA and a pitch of 1mm. The captured 3D images have slice thickness of 0.75mm with typical pixel spacing of 0.6x0.6 mm and slice separation of 0.5 mm.

Bronchial geometry acquisition
The airway analysis was performed with the Pulmonary Workstation 2 (PW2) software (VIDA Diagnostics, Coralville, IA). The lungs, lobes and airway tree bronchi were segmented from the HRCT images down to the most distal distinguishable generation. A set of measurements was then performed on each segmented bronchus, including the average Din, Ai, WT, Pi and Po (outer perimeter). The manual and semi-automatic use of PW2 as a research tool is explained in detail in [36]. The airways in each segmented image were then uniquely identified, anatomically labeled and matched (Fig 1); only bronchi that were available at both levels of respiration were used. Standard nomenclature was used, matching the anatomical airway names with their respective acronyms [52]. Briefly, the airway names and their acronyms are the following: T (trachea), R (right main bronchus), L (left main bronchus), RMI (right intermediate), RU (right upper lobe), RM (right middle lobe), LU (left upper lobe) and LL (left lower lobe). The right lung segmental bronchi are for the right upper lobe RB1 (apical), RB2 (posterior), RB3 (anterior), for the right middle lobe RB4 (lateral) and RB5 (medial) and for the right lower lobe RB6 (superior), RB7 (medial basal), RB8 (anterior basal), RB9 (lateral basal), RB10 (posterior basal). For the left lung, the segmental bronchi in the upper lobe are LB1 (apical), LB2 (posterior), LB3 (anterior), LB4 (inferior), LB5 (superior) and in the lower lobe LB6 (superior), LB8 (anterior medial basal), LB9 (lateral basal) and LB10 (posterior basal).
The Dh was calculated according to the formula while airway distensibility was defined as where Ai TLC and Ai MLV is the airway lumen area at TLC and MLV respectively and V TLC and V MLV are the volumes of the bronchi's subtended lungs (right or left depending on bronchus position). Paired Student t-tests were used to examine statistical significance between parameters (Dh vs Din) at both lung volumes (p<0.05). Linear regression analyses where assessed using the Pearson coefficient. For practical purposes all lobar bronchi were assigned generation 2, segmental bronchi generation 3 etc. Given this definition, airway segmentation of matched airways was reliably available down to generation 5 only in the RB1, RB3 (right upper lobe) RB6, RB10 (right lower lobe) LB3 (left upper lobe) and LB10 (left lower lobe) segments.

Results
The results from the Dh calculations are provided in Table 2 for both MLV and TLC groups. For convenience these results were organized per generation and per lobe. All bronchi displayed a 13.5±9% positive increase in Dh on average, ranging between 8.5±4.7% in the trachea and main bronchi and 19.6±9% in generations 4-5 of the left lower lobe. The lower lobes displayed higher % change on average than the upper lobes. In generation 3 that difference is 4.75% and 3.1% for the right and left lung respectively, while this trend persists in generations 4 (4.9% and 6.8%) but not fully in generation 5 (-0.7% and 6.2%). A decrease in Dh was also observed after full inspiration, but was rare (<5% of the examined samples). A Student t-test between the two groups revealed that, with the exception of RL lobe generations 2 and 5, all other comparisons were statistically significant (p<0.05).
A linear regression analysis of Dh between MLV and TLC shows that the majority of the lobar regions display moderate or strong (R 2 >0.55) correlation to lung volume change (see Fig  2). In fact, the lowest correlation coefficients can be observed at generation 5 for the RU, RL and LU lobes, where measurement artifacts and more sparsely populated datasets are common due to diminishing airway sizes. The graphs of Fig 2 represent the effect of the two lung volumes on the Dh for the trachea, right and left main bronchi (Fig 2A), the lobar bronchi ( Fig  2B), and the bronchi of generations 3-5 for the RU lobe ( Fig 2C, 2D and 2E). With very few exceptions, the Dh increases for the large majority of the samples involved, with the trend lines almost parallel to the identity line, thus verifying the 13.5% average increase in Dh between lung volumes. These results were similar for all lobes.
The comparison between Dh and Din, displayed in Fig 3, gives statistically significant results for all the airways involved in this study, with Din>Dh and very high correlation coefficients, as the R 2 of the linear regression analyses are between 0.8734 ( Fig 3B) and 0.9683 ( Fig  3C). Very similar results could be observed for the other lobes for both MLV and TLC. The average ratio Din/Dh was calculated to be 1.124±0.068 at MLV and 1.108±0.063 at TLC. Even though only 9 out of the 52 examined bronchi displayed statistically significant differences between the two lung volumes (significance given by p<0.01), almost all individual bronchial segments showed that the ratio Din/Dh is larger in MLV than in TLC.
The Ai and Pi were collected in Table 3 for all uniquely identifiable segmented bronchi in MLV and TLC (Ai in mm 2 and Pi in mm). It can be seen that only 6 out of 19 possible generation 5 bronchial regions could be analyzed in the MLV images. For that reason, comparisons between the two lung volumes are only possible in these 5 particular regions. The lower lobes display higher change in airway caliber after inflation to TLC when compared to the upper lobes. For generation 3, the Ai change was 3.3% and 4.2% for the right and left lungs, respectively, while similar trends were found for generations 4 (6% and 13.4%) and 5 (4.9% and 12.2%). The average distensibility per airway, as defined in Eq (2), is recorded in Table 4. A box and whisker plot of these data for all available bronchi of generations 3 and 4 can be seen in the supplemental graph of S1 The correlation between the lumen area change and lung volume change, as expressed in the same equation, was also examined and recorded in Table 4 in the form of the correlation coefficient R 2 . The plotted results of these correlations can be seen in   ( Fig 4A). Equivalent trends could be observed for the same lung region for the generation 4 data (Fig 4B) with somewhat higher R 2 values (R 2 RB2-gen4 = 0.4444). Fig 4C and 4D show generation 3 and 4 data respectively for the right lower lobe. Airways RB8-10 show trends similar to the ones observed in the right upper lobe. An examination of generation 4 data indicates that, with the exception of airway RB10 (R 2 RB10-gen3 = 0.589), the change of lung volume does not correlate very well with the changes in lumen area. It can be observed that the distensibility of bronchus RB7 is consistently lower than that of RB6, 8, 9 and 10, even though statistical significance is only reached at generation 4 (p<0.05). The excessive variability of airway RB6 is also notable.
A lower dispersion of data can be observed in the left lung for all airways. Many analyzed bronchi, including LB3 (gen. 3 and 5), LB4, LB5, LB6 (gen. 4) LB9 (gen. 3 and 4) and LB10 (gen. 3, 4 and 5) had R 2 > 0.4. Overall, the left lung had higher R 2 values when equivalent regions from both lungs were compared. The influence of the left lung volume is distinct but relatively similar for bronchi LB8, LB9 and LB10 (Fig 4E). The generation 4 data can be seen in Fig 4F. Furthermore, statistical analysis revealed greater intra-region variability as compared to the right lung. For example, airway LB3 has lower distensibility (p<0.05) than both LB1 and LB2 in generation 3, while most comparisons between generation 4 bronchial distensibilities also gave statistically significant differences. As opposed to RB6, LB6 has lower variability and has the highest distensibility of all generation 3 airways of the left lung (p<0.05). Finally, similar levels of distensibility are observed between generations 3 and 4 in the right lung, while in the left lung statistical significance is reached only for airways LB1 (gen3>gen4, p = 0.0041) and LB3 (gen3<gen4, p = 0.0419). The Pearson correlation coefficient was found to be R 2 = 0.4793 (full dataset) and R 2 = 0.7401 (LB1 excluded). The distensibility of generation 3 vs generation 4 can be seen in Fig 5.

Discussion
Even though several studies in the past have investigated the morphometry of the human bronchi during deep inspiration, only one study was found containing data on the Dh, calculated for airways at TLC and normalized to predicted trachea diameter [38]. If the TLC data of Table 2 are normalized to the value of the trachea for the purpose of comparison with Choi et al. [38], it is observed that our Dh values are within their overall range covering healthy, non severe and severe asthmatic patients, but do not strictly match the asthmatic range (S1 Table). Airway shape and geometry in different lung volumes in asthma Also, somewhat larger discrepancies appear in the left lung. This could be attributed to the known measurement artifacts in the left lung due to heart motion. Another explanation could be that the results of the Mch challenge administered to the patients in our study induce higher spatial variability within the lung in terms of broncho-constrictive effects.
Our results indicate that lower lobe bronchi change more in caliber than the upper lobe airways (Tables 2 and 3, Fig 4) at maximal lung inflation, especially in the left lung. In general, differences between lung regions have been reported in the literature for both healthy subjects [53] and COPD patients [54]. In the case of asthma, this difference may indicate that the effects of the methacholine challenge are more prominent in the upper regions of the lung. For example, Kotaru et al. [55], in a study examining the possible locations and patterns of airway response after different types of broncho-provocation on asthmatic subjects, did not report any regional differences but did find significant variability, especially with respect to the airway depth where the challenge is more effective. Furthermore, they observed complex morphometric alterations, not uniform nor graded, with mixed pattern dilation-constriction within and among anatomical regions. Another possibility is that the upper lobes of the human lung are more prone to disease, similar to COPD, as opposed to the lower lobes, the latter exhibiting more endurance during challenge but being faster to decline in function when more serious symptoms are observed [56]. This would be in line with the findings of [38], i.e. preferential regional occurrences of wall thickening in the upper lobes and airway narrowing in the lower lobes.
In addition to the above, larger diameter and luminal area changes between MLV and TLC exist in the more distal airways of the left lung (generations 4-5) rather than more proximally (generations 2-3, Tables 2 and 3). While investigating the effects of inspiration on airway dimensions of healthy subjects, Petersen et al. [53] pointed out that airways distend fully from generation 5 and deeper since their tissue is dominated by smooth muscle. Hoshino et al. [37] confirmed that there is an increasing correlation between Ai and clinical indices of disease as we penetrate deeper into the lung, confirming that airway remodeling takes place more distally. A similar study by Shimizu et al. [57] had comparable conclusions. Our results seem to Airway shape and geometry in different lung volumes in asthma partially agree with those of Brown et al. [58], who detected higher distensibilities (measured as %Ai change) in airways with diameter <3mm after Mch challenge. This could be attributed to the higher edema present in the smooth muscle, which may attenuate the radial traction produced by the inflation [59] or to airway remodeling of the extra-cellular matrix [45].
The high correlation coefficients and combined data regression lines of Fig 2 demonstrate a relatively mild and stable linear relationship for Dh between MLV and TLC, persistent throughout the lungs with some minor local variations. It is interesting to note that there is a relatively small number of bronchi displaying either no change or even a reduction of Dh at TLC (samples below the identity line). This dilation failure, or even the constriction, after inspiration at TLC might be statistically insignificant appositely to airway tree morphology but being located at proximal bronchi can have a considerable effect on symptoms [55,60]. Recently, Kim et al. [10] while investigating a coupled model of tissue deformation and network airflow observed that small airway constrictions can disturb flow patterns over large areas, reinforcing the findings of Venegas et al. [28,61] who showed that sudden closures of small airways break the system balance, leading to catastrophic shifts. They also found that the constriction in large airways that feed large volumes can profoundly affect lung function, an event that is not uncommon as we have discovered in this and a previous study by our group on asthmatic patients [36], where severe broncho constrictions were discovered as proximally as generation 3. The ventilation heterogeneity observed in asthma is therefore connected to regional anatomic heterogeneity, a fact confirmed by [7], making the use of average values over an entire generation misleading. Unfortunately, no method for the prediction of location and degree of proximal airway constrictions is currently available. Thus, morphometric models of disease that fail to give accurate representations of these defects, both geometrically and spatially, would result in functional models that could not exhibit important phenomena.
The comparison between the Dh and the Din demonstrated that the latter is almost uniformly larger than the former (Fig 3) by 12.5% at MLV and 10.8% at TLC, with the difference between lung volumes displaying statistical significance (p<0.01). In some airways, this difference can be considerably higher; for example, it was found to be 23±18% in the RB6 generation 3 airway. While these differences in diameter are relatively minor, functionally it should be remembered that for laminar and fully developed flow in a tube, the pressure drop is directly proportional to the flow rate but inversely proportional to the fourth power of diameter. Thus, a decrease in diameter by 10% would result in a 52% increase in pressure drop, all else being equal [48]. It is unclear if this systematic difference is a result of the measurement method used by the segmentation software or if it reflects a physical reality of the diseased human lung. If non-circular sections are considered throughout the segmented airway tree, as is implied by our results and also by [18,38], finding the Dh should be a priority over the measurement of Din for reasons concerning both airway shape and overall airflow resistance.
Since airway resistance depends on the Dh [48], the effect of this difference has a multiplicative effect on overall tree resistance calculations as the number of airways doubles with each successive generation distally into the lung [62]. Furthermore, even though the Din/Dh displayed statistically significant differences between MLV and TLC in only 9 out of 52 examined segments, the ratio was found to be consistently larger at MLV, indicating a trend for higher airway resistance at this breathing condition. The analysis of HRCT, due to the limited spatial resolution, was not sufficient to reveal the full complexity of the bronchial perimeter shape, providing an elliptical cross-section approximation. For that reason, the difference in the Din/ Dh ratio between the two lung volume states probably represents a small increase in airway circularity at TLC; i.e., there is an anisotropic airway expansion across its circumference during inspiration. This is probably related to the elastic properties of the lung parenchyma, which is an important feature of the mechanics of breathing in disease conditions [10,63].
Several previous studies have reported the bronchial lumen area, either directly or as percentage of total area, for various uniquely identified airways and for various states of health and disease. The results were mostly presented either as average values per airway generation [38,64] or are given only for specific, easily identifiable generation 3 airways that may act as surrogate markers of lung function and tree morphology [34,54,[65][66][67], or both [37,53,57]. In this study, Ai and Pi measurements were provided in Table 3 for all segmented bronchi for generations 3-5 for both MLV and TLC and without any normalization. This was done mostly because, as mentioned earlier, concentrating on specific airways may ignore localized characteristics of disease or disease progression due to the intense intra-subject variability of asthma. Furthermore, as Bakker et al. [54] suggested, bronchial measurements should not be presented as an average value for the whole lung, a statement supported by the findings of Fain and coworkers [68] who found that normalization of airway measures did not eliminate dependence on airway segment and other patient specific factors in a group of normal and asthmatic subjects. Finally, individualized airway morphometric information can be really useful for lung modeling applications [11,21] Many techniques have been used in the past to determine what is the influence of lung inflation, if any, on the airway geometry [13,54,58,69]. In this study, a definition of airway distensibility that explores relative changes in Ai and ΔV was adopted. The results are reported for all bronchial segments in Table 4. The influence of volume change per lung on Ai change was also reported as the Pearson correlation coefficient and is displayed in the graphs of Fig 4. These data suggest that volume changes in the right lung only moderately affect the Ai in the relevant airways, while better correlations exist when the left lung is considered. This might be due to the fact that the morphometric measurements conducted in this study where made on static images (HRCT) while the airway distensibility has been shown to be primarily affected by the dynamic ability to change lung volumes in asthma [70].
The formulas of distensibility suggested by Bakker et al. [54] and Diaz et al. [35] were also examined. Even though the best correlations between ΔAi and ΔV were observed with the definition adopted by [35], we believe that the definition of Eq (2) makes the assessment of bronchial distensibility easier both within and between different airways and different generations. This is because, if the value of distensibility in Eq (2) is unity, one can claim that the change in area is the same as the change in lung parenchyma, i.e. the airway and the parenchyma display similar compliance. In the group of patients analyzed in this study, the distensibility was calculated almost unilaterally below unity, with an average of 0.45±0.28, displaying variation between different lung regions. This failure of the central airway lumen to dilate to an equivalent degree with the inflation of the lung should be expected due to the tissue composition of these bronchi. For example, [53] considers that full airway distensibility is only achievable after generation 5 in healthy lungs, where soft tissue is dominant. However, in asthma extra bronchial rigidity might be expected due to mechanical influence of the bronchial wall edema,. Furthermore, airway wall contractility should also be affected by the increased tone effected by the methacholine challenge. Unfortunately, no healthy subjects were included in this study for comparison purposes. It has to be noted that if the effects of respiration and lung volumechange on airway size are to be comparable across studies, similar definitions must be universally adopted.
In the context of mathematical modeling of the bronchial tree, the results provided in this paper can potentially be valuable, especially when parent-daughter relationships are used to calculate bronchial geometry information of distal lung regions [2,8,11,41]. The high correlation between Din and Dh allows for a straightforward linear conversion while the data of Tables 2, 3 and 4 provide sufficient information on Dh, Ai and Pi for the relevant adjustments to be made to the models, depending on respiratory conditions, initial information and other system requirements for the modeling application. Many modern tracheo-bronchial tree models utilize medical images of the lung acquired at TLC to obtain the initial conditions for airway tree propagation to the acinar level [11,21,71]. These models, in turn, are used to simulate conditions at regular breathing without the geometric (circularity, resistance, non-persistent stenoses) and ventilation discrepancies being taken into account. This misconception has recently been recognized by researchers using such models of the lung. For example, in a CFD study on aerosol deposition in bronchial tree models extracted from HRCT images at MLV and TLC, Katz et al. [27] concluded that lung volume should be explicitly considered when using deposition models. Bordas et al. [21], lacking relevant information for a population of asthmatic patients, used the minimum inner diameters of the segmented airways in conjunction with a circular cross-section assumption to predict airflow in the lung. Another method was used by Kim et al. [10], who implemented their dynamic ventilation model only after uniformly scaling down the bronchial tree model obtained at TLC by a scale factor ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi FRC=TLC 3 p . One limitation of this study is the lack of sufficient resolution in the medical imaging modality used to assess the real shape of airway internal cross-section perimeter. Especially in disease, this perimeter might have a complex geometry due to swelling, hyperplasia of goblet cells, mucous excretions etc. [72] which, in turn, could considerably affect airflow and particle deposition. Therefore, these changes in circularity between MLV and TLC that are reported here could be connected to the "unfolding" of the bronchus as the pressure changes in the lung drive the airways to obtain higher volumes. Furthermore, comparisons between MLV and TLC could only concentrate on airways that could be segmented in both lung volumes. This undoubtedly biases the data in favor of larger airways, especially in more distal lung regions where image resolution is not sufficient for airway segmentation at MLV. Finally, as stated earlier, no dynamic effects of respiration could be taken into account as there is no temporal component in the medical images used herein.
In conclusion, this study was concerned with the examination of the hydraulic diameter human bronchial tree as calculated from the lumen area and inner perimeter at both mean lung volume and total lung capacity. The relationships between several morphometric features were investigated, along with the effects of lung volume change on airway shape. The airway hydraulic diameter was found to be consistently higher than average inner diameter at both investigated lung volumes, a fact that can have considerable effects on resistance and pressure loss at proximal lung regions. Furthermore, differences in caliber changes were observed between upper and lower lobe bronchi, which might be indicative of differences in asthma expression in different lung regions. A failure of some airways to dilate at TLC was also documented, confirming the mixed constriction-dilation patterns of asthmatic airways that lead to the patchiness of lung ventilation characteristic of this disease. Finally, generation 3 distensibility was predictive of distensibility at generation 4. These results could be used for adapting models of the human trachea-bronchial tree and the interpretation of results for ventilation and aerosol deposition simulations. The probability plot comparing the distensibility data to the normal distribution (reference line--), created for the right and left lung generation 3 and 4 data. It can be observed that the fit of the data to the reference line deteriorates when distensibility values exceed one. (TIF) S1 Table. Normalized hydraulic diameter data. The hydraulic diameter data are normalized to the hydraulic diameter calculated for the trachea and are collected for the same airway tree bronchial regions as the ones considered by Choi et al. [38] (where available). The data are give as average values ± standard deviation values. (DOCX) S1 File. Data processing file. This is the excel file used to gather and statistically analyze the bronchial tree morphometric information that are investigated in this study. The file contains the average inner diameter, average lumen area, average inner perimeter, calculated hydraulic diameter and several definitions for airway distensibility. The data are gathered in a variety of ways including per individually identified airway and per generation. The file also contains statistical analysis.