1. Motivation
Asthma and chronic bronchitis are the most prominent obstructive lung diseases affecting millions of people worldwide. Asthma affects 8% of the adult population and 20% of children [
1]. Chronic bronchitis, a manifestation of chronic obstructive pulmonary disease, affects 25% of the adult population and is the fourth leading cause of death [
2]. Asthma is generally classified as
reversible, triggered by allergens, weather or exercise [
3], while chronic bronchitis is
irreversible, triggered by pollutants, toxins and smoke [
4,
5].
Figure 1 illustrates the two common features of asthma and chronic bronchitis:
airway constriction associated with smooth muscle thickening at the outer airway wall; and
airway inflammation associated with mucosal growth at the inner airway wall [
6]. These two phenomena manifest themselves mechanically in an increase in pressure from active smooth muscle cell contraction and an increase in volume from an influx of mononuclear cells flooding the inner wall lining. Collectively, both phenomena initiate an inward folding of the airway wall, ultimately resulting in a narrowing of the lumen and progressive airflow obstruction [
7,
8].
Figure 1.
Airway wall remodeling in chronic lung disease. In contrast to the healthy airway wall, the airway wall in asthma, left, and chronic bronchitis, right, displays
airway constriction, associated with smooth muscle thickening at the outer airway wall, and
airway inflammation, associated with mucosal growth at the inner airway wall. These changes manifest themselves in an increase in pressure at the outer wall and an increase in volume at the inner wall, resulting in a narrowing of the lumen and progressive airflow obstruction; adapted from [
1,
4].
Figure 1.
Airway wall remodeling in chronic lung disease. In contrast to the healthy airway wall, the airway wall in asthma, left, and chronic bronchitis, right, displays
airway constriction, associated with smooth muscle thickening at the outer airway wall, and
airway inflammation, associated with mucosal growth at the inner airway wall. These changes manifest themselves in an increase in pressure at the outer wall and an increase in volume at the inner wall, resulting in a narrowing of the lumen and progressive airflow obstruction; adapted from [
1,
4].
Figure 2 displays a schematic of the human lungs. The trachea, the root of the respiratory tree, branches into two main bronchi, which enter the left and right lungs. In the lungs, the bronchial tree continues to branch in humans into 23 to 27 generations, yielding approximately 17 million branches [
9]. The final generation of terminal bronchioles opens into the alveolar space, the region where gas transfer takes place [
6]. Airway wall remodeling affects the small airways between the fourth and fourteenth generation [
4]. Healthy small airways are less than 2 mm in diameter, non-cartilaginous and compliant [
10]. In chronic lung disease, the dimensions of the small airways roughly double, resulting in a loss of structural compliance, airflow obstruction and difficulties in breathing [
11].
Figure 2.
Schematic of the lungs. The trachea, the root of the respiratory tree, branches into two main bronchi, which enter the left and right lungs. In the lungs, the bronchial tree continues to branch in humans into 23 to 27 generations, yielding approximately 17 million branches. The final generation of terminal bronchioles opens into the alveolar space, where gas transfer occurs. Airway wall remodeling affects the small airways between the fourth and 14th generation. Healthy small airways are less than 2 mm in diameter, non-cartilaginous and compliant. In chronic lung disease, the small airways experience smooth muscle thickening and mucosal growth, resulting in progressive airflow obstruction.
Figure 2.
Schematic of the lungs. The trachea, the root of the respiratory tree, branches into two main bronchi, which enter the left and right lungs. In the lungs, the bronchial tree continues to branch in humans into 23 to 27 generations, yielding approximately 17 million branches. The final generation of terminal bronchioles opens into the alveolar space, where gas transfer occurs. Airway wall remodeling affects the small airways between the fourth and 14th generation. Healthy small airways are less than 2 mm in diameter, non-cartilaginous and compliant. In chronic lung disease, the small airways experience smooth muscle thickening and mucosal growth, resulting in progressive airflow obstruction.
Figure 3 illustrates a typical cross-section of the small airway wall with its three distinct layers. The mucosa, the innermost layer, consists of the epithelium, the basement membrane and a subepithelial collagen layer. The stiffness of the mucosa depends primarily on the subepithelial collagen layer, whose thickness ranges from 4–5 μm in healthy to 7–23 μm in diseased airways [
12]. The submucosa, the middle layer, consists of fibroblasts and proteoglycans embedded in a loose, irregular network of elastin and collagen. Although the submucosa is about an order of magnitude softer than the mucosa [
6], it contributes significantly to the overall stiffness of the airway wall because of its considerable thickness. The smooth muscle layer, the outermost layer, consists of spirally arranged smooth muscle cells. In the healthy airway, these smooth muscle cells regulate the amount of airflow through active cellular contraction. In the diseased airway, smooth muscle thickening and uncontrolled muscle contraction generate an elevated contractile force, resulting in an increased wall stiffness and a decreased lumen.
Figure 3.
Schematic of the airway wall. Small airways consist of three distinct layers, the mucosa, the submucosa and the smooth muscle layer. In chronic lung disease, the smooth muscle layer thickens and creates an elevated pressure at the outer wall, while the mucosal layer experiences inflammation and increases in volume at the inner wall. The airway wall folds inwards, resulting in progressive airflow obstruction.
Figure 3.
Schematic of the airway wall. Small airways consist of three distinct layers, the mucosa, the submucosa and the smooth muscle layer. In chronic lung disease, the smooth muscle layer thickens and creates an elevated pressure at the outer wall, while the mucosal layer experiences inflammation and increases in volume at the inner wall. The airway wall folds inwards, resulting in progressive airflow obstruction.
The numerous fatalities and sheer number of people whose lives are impacted daily by asthma and chronic bronchitis motivate the need to understand the microstructural mechanisms of chronic lung disease. The past two decades have seen a tremendous progress in the mathematical and mechanical analysis of airway constriction and airway inflammation [
6,
13].
Table 1 classifies the most prominent analytical and computational models of chronic lung disease. While earlier approaches focus primarily on airway constriction and smooth muscle thickening modeled through an increased pressure of the outer layer [
14,
15], recent approaches focus on airway inflammation and mucosal thickening modeled through an increased volume of the inner layer [
8,
16]. Few studies combine both mechanisms to explore whether the associated structural alterations have a positive or negative feedback on one another [
13,
17].
Analytical modeling can provide valuable insight into the critical failure mode [
18,
19]. Using classical bifurcation analysis [
20], analytical models can help identify critical conditions, such as the critical pressure [
14,
21] or the critical amount of growth [
8,
16] at the onset of failure. However, analytical models fall short of predicting the progression of failure throughout the post-failure regime [
22,
23]. This is particularly important, since complications in chronic lung disease are typically associated with the later stages of failure rather than with the onset of folding. In addition, analytical models are typically limited to simplified constitutive models and regular geometries. It is thus not surprising that all existing models of airway wall thickening are restricted to regular, circular cross-sections [
24].
Numerical modeling cannot only predict the onset of failure, but also the structural response throughout the entire post-failure regime [
25,
26]. For example, they can predict the critical growth at first contact, which serves as a valuable clinical metric to characterize the severe stages of airway obstruction. Using finite element analysis, numerical models can easily incorporate more realistic constitutive models and complex patient-specific geometries [
27,
28]. Computational modeling has been intensely used to simulate growth within the vascular tree [
29,
30]. Yet, to date, there are no compelling computational models to simulate growth within the bronchial tree. This is particularly important, since the bronchial tree branches a lot more frequently than the vascular tree. In fact, airway branching occurs once every two to four airway diameters [
6]. This implies that only a small fraction of the bronchial tree is in fact a circular cross-section, and irregular, non-circular cross-sections might play an important, yet underestimated, role.
Table 1.
Classification of existing models of chronic lung disease.
Table 1.
Classification of existing models of chronic lung disease.
analysis | microstructural mechanism | analytical solution onset of failure | numerical solution post-failure regime |
---|
symptoms |
---|
airway constriction | smooth muscle thickening → increased pressure on submucosal layer | bifurcation analysis → folding pressure → number of folds [13,14,15,17] | finite element analysis → contact pressure → folded configuration [14,15] |
airway inflammation | mucosal inflammation → volume growth of mucosal layer | bifurcation analysis → folding growth → number of folds [8,13,16,17,23,24] | finite element analysis → contact growth → folded configuration [8,16,23,25,31,32] |
We hypothesize that the geometry of the bronchial tree plays a crucial role in chronic airway obstruction and that critical failure conditions vary significantly along a branching airway segment. To test this hypothesis, we introduce a continuum model for mucosal thickening based on the theory of finite growth [
33,
34] and establish a computational model for its numerical solution using nonlinear finite element analysis. We represent growth through the multiplicative decomposition of the deformation gradient [
35] into an elastic tensor, which induces stress [
36], and a growth tensor, which we prescribe constitutively [
37,
38]. To mimic chronic disease progression up to and beyond the onset of airway wall folding, we gradually increase the amount of growth [
31,
39]. Although the underlying concept is generally applicable to both isotropic and anisotropic growth [
40,
41], here, we assume that growth is purely isotropic. We represent it through a spherical tensor scaled by the amount of volume growth [
42]. We embed the resulting growth model into a nonlinear finite element setting and solve it using an incremental iterative Newton–Raphson solution scheme [
43,
44]. In the sequel, we outline the details of our model.
The remainder of this manuscript is organized as follows: We briefly summarize the continuum and computational modeling of finite growth in
Section 2. To identify critical regions of airway narrowing, we perform systematic studies along a branching airway segment in
Section 3: First, we visualize the rich variety of possible failure modes along a branching airway segment by inducing prescribed failure modes in elliptical cross-sections through sinusoidal perturbations in
Section 3.1. Next, we explore the role of the relative mucosal thickness by inducing natural failure modes in circular cross-sections through random perturbations in
Section 3.2. Last, we explore the role of the relative mucosal thickness by inducing natural failure modes in elliptical cross-sections without perturbations in
Section 3.3. We compare the results against the literature and discuss the limitations of our study in
Section 4. Finally, we conclude with a brief discussion in
Section 5.
3. Results
To identify critical regions of airway narrowing, we perform systematic parameter studies for varying cross-sections within the bronchial tree.
Figure 4 illustrates four representative elliptical cross-sections along a branching airway segment. In analogy with the literature [
8,
15], we introduce the airway radius
R as the distance from the airway center to the closest point on the mucosal-submucosal interface. We denote the mucosal and submucosal thicknesses as
and
, such that the inner and outer radius of a circular cross-section are
and
. For elliptical cross-sections, we characterize the degree of ellipticity through the semi-major and semi-minor axes
and
of the submucosal-muscular interface, where the semi-minor axis is identical to the outer radius of submucosal layer
. We model the mucosal and the submucosal layers as elastically incompressible with stiffnesses
and
and the idealize the smooth muscle layer as rigid. While the mucosa is allowed to grow at a speed of
τ towards a maximum growth of
, the submucosa remains purely elastic.
Figure 4.
Branching airway segment. Along the bronchial tree, the airway cross varies substantially in geometry. We characterize each cross-section through the semi-major-to-semi-minor-axis ratio and through the mucosal-thickness-to-radius ratio . The ratio is larger close to and smaller away from a branching region. The ratio is smaller close to and larger away from the trachea.
Figure 4.
Branching airway segment. Along the bronchial tree, the airway cross varies substantially in geometry. We characterize each cross-section through the semi-major-to-semi-minor-axis ratio and through the mucosal-thickness-to-radius ratio . The ratio is larger close to and smaller away from a branching region. The ratio is smaller close to and larger away from the trachea.
To characterize the individual cross-sections in terms of non-dimensional parameters, we introduce the mucosal-thickness-to-radius ratio
, the submucosal-thickness-to-radius ratio
, the ellipticity ratio
, and the mucosal-to-submucosal stiffness ratio
. Since previous studies have shown that the failure mode is more sensitive to the ratio
than to the ratios
and
[
8,
15], we fix the latter two to
and
[
6]. In the following studies, we vary the relative mucosal thickness ratio between
and the ellipticity ratio between
.
We discretize each cross-section using classical tri-linear brick elements and assume plane strain conditions along the airway’s long axis. Since the smooth muscle layer is significantly stiffer than the mucosal and submucosal layers, we represent it implicitly through homogeneous Dirichlet boundary conditions at the submucosal-mucosal interface. We allow the mucosal layer to grow and fold inward. Once we detect the first contact between two folds, we terminate the simulation and record the corresponding growth factor as the critical growth at first contact. To explore regional variations of growth along the bronchial tree, we perform three sets of simulations: In
Section 3.1, for varying elliptical cross-sections with
sinusoidal perturbations, we demonstrate how growth drives the solution along the prescribed failure mode, which is
a priori defined by the corresponding perturbation. In
Section 3.2, for circular cross-sections and
random perturbations, we demonstrate how growth drives the solution into the natural failure mode, which is governed by the relative mucosal thickness. In
Section 3.3, for varying elliptical cross-sections
without perturbations, we demonstrate how growth drives the solution into the natural failure mode, which is induced naturally by curvature heterogeneity and governed by the relative mucosal thickness.
3.1. Sensitivity of Failure Mode with Respect to Ellipticity
The first set of examples in
Figure 5 visualizes the wealth of failure modes for different ellipticity ratios. From a mechanical point of view, this study is closely related to the buckling of elliptical tubes under radial compression [
45]. We explore six different cross-sections with ellipticity ratios varying as
, where the limits of
and
mimic cross-sections away from and close to a branching region. To trigger a specific,
prescribed failure mode, we perturb the homogeneous elliptical cross-section with
sinusoidal perturbations of a magnitude of one fourth of the mucosal thickness, 0.25
, and a frequency of
modes along the inner mucosal surface. To visualize potential failure modes, we vary the mode number as
and perform two simulations for each mode: Simulation I is at least symmetric to the semi-major axis and has at least one fold along the axis of symmetry; Simulation II is at least symmetric to the semi-minor axis.
Figure 5.
Sensitivity of the failure mode with respect to ellipticity. Smaller ellipticity ratios represent circular cross-sections away from a branching region; larger ellipticity ratios represent elliptical cross-sections close to a branching region. With increasing ellipticity, the heterogeneity of the failure mode increases. With increasing heterogeneity, the lumen area at the first contact remains significantly larger. This might indicate that circular cross-sections away from the branching region are at a higher risk of to airflow obstruction than elliptical cross-sections close to a branching region.
Figure 5.
Sensitivity of the failure mode with respect to ellipticity. Smaller ellipticity ratios represent circular cross-sections away from a branching region; larger ellipticity ratios represent elliptical cross-sections close to a branching region. With increasing ellipticity, the heterogeneity of the failure mode increases. With increasing heterogeneity, the lumen area at the first contact remains significantly larger. This might indicate that circular cross-sections away from the branching region are at a higher risk of to airflow obstruction than elliptical cross-sections close to a branching region.
Figure 5 illustrates the different failure modes for varying ellipticity ratios
and varying mode numbers
. The color code reflects the circumferential stress with red values at the inner mucosa corresponding to high circumferential traction and blue values at the outer submucosa corresponding to low traction. As the ellipticity ratio increases, from left to right, the heterogeneity of the failure mode increases. Notably, all folding patterns, from top to bottom, display a discontinuous alteration in the underlying failure mode as the ellipticity increases: Mode
, displayed in the first row, changes from a double horizontal double vertical contact mode at
via a double horizontal contact mode at
to a double vertical contact mode at
. Mode
, displayed in the second row, changes from a quadruple-diagonal contact mode at
via a single vertical contact mode at
and
to a double vertical contact mode at
. Mode
, displayed in the sixth row, changes from a sextuple contact mode at
via a double horizontal contact mode from
to
to a double vertical contact mode at
. Mode
, displayed in the eighth row, changes from a septuple contact mode at
via a single vertical contact mode from
to
to a double horizontal contact mode at
. With increasing ellipticity, from left to right, the lumen area at the first contact remains significantly larger. This might indicate that circular cross-sections away from the branching region are at a higher risk of airflow obstruction than elliptical cross-sections close to a branching region.
3.2. Sensitivity of the Failure Mode with Respect to Relative Mucosal Thickness
The second set of examples in
Figure 6 probes the sensitivity of the failure mode for different relative mucosal thickness ratios in a circular cross-section with
. From a mechanical point of view, this study is closely related to the buckling analysis of circular tubes under an external pressure [
17]. We explore eight different cross-sections with relative mucosal thickness ratios varying as
. As illustrated in
Figure 4, variations of
occur not only in similar generations of bronchi with the same radius
R and varying initial mucosal thickness
, but also in different airway generations along the bronchial tree with varying radius
R and similar mucosal thickness
. In fact, the overall wall-thickness-to-radius ratio was measured to vary from 0.044 in the first generation of bronchi to 0.075 in the sixteenth generation [
6], indicating that the relative thickness increases as we descend down the bronchial tree. To trigger the
natural failure mode, we perturb the homogeneous circular cross-section with
random perturbations of one percent of the mucosal thickness, 0.01
. To mimic the onset of failure, we take snap shots at
, which correspond to the first points of contact of the eight different cross-sections.
Figure 6 illustrates the evolution of the different failure modes for varying relative mucosal thicknesses
and for different first points of contact
ϑ. The upper right section displays the circumferential stress; the lower left section displays the three-dimensional failure mode, both synchronized at the same time point.
Figure 6 demonstrates that the relative mucosal thickness determines the failure mode and the number of folds, irrespective of the small initial random perturbation. As growth progresses, the folds become more and more pronounced until first contact occurs.
Figure 6.
Sensitivity of failure mode with respect to relative mucosal thickness for a circular cross-section with . Larger relative mucosal thicknesses represents smaller bronchi away from the trachea; smaller relative mucosal thicknesses represents larger bronchi close to the trachea. With increasing relative mucosal thickness, the number of folds decreases. With the decreasing number of folds, the critical growth ϑ at the first contact becomes significantly smaller. This might indicate that smaller bronchi with a larger relative mucosal thickness are at a higher risk of airflow obstruction than larger bronchi with a smaller relative mucosal thickness.
Figure 6.
Sensitivity of failure mode with respect to relative mucosal thickness for a circular cross-section with . Larger relative mucosal thicknesses represents smaller bronchi away from the trachea; smaller relative mucosal thicknesses represents larger bronchi close to the trachea. With increasing relative mucosal thickness, the number of folds decreases. With the decreasing number of folds, the critical growth ϑ at the first contact becomes significantly smaller. This might indicate that smaller bronchi with a larger relative mucosal thickness are at a higher risk of airflow obstruction than larger bronchi with a smaller relative mucosal thickness.
Figure 6 confirms that circular cross-sections with random perturbations can generate both even and odd failure modes. As the relative mucosal thickness decreases, the number of folds increases from
at
to
at
. For cross-sections with relatively thin mucosal layers, towards the right column and the bottom row, more growth is required to form the first contact. The overall mucosal thickness at this first contact point, however, is remarkably similar for all eight cross-sections, as shown along the diagonal. In summary, with increasing relative mucosal thickness, from top right to left, the required growth to form the first contact decreases. This might indicate that smaller bronchi with larger relative mucosal thickness, displayed towards the left, are at a higher risk of airflow obstruction than larger bronchi with a smaller relative mucosal thickness, displayed towards the right.
3.3. Sensitivity of Failure Mode with Respect to Ellipticity and Relative Mucosal Thickness
The third set of examples in
Figure 7 and
Figure 8 probes the sensitivity of the failure mode for both different ellipticity ratios and different relative mucosal thickness ratios. We compare a moderately elliptical cross-section with
and a severely elliptical cross-section with
to the circular cross-section with
discussed in the previous subsection. For each cross-section, we vary the relative mucosal thickness as
. Elliptical cross-sections display a natural heterogeneity as the curvature varies along the circumference. To trigger the
natural failure mode, we can thus simulate the plain elliptical cross-section
without perturbations. To mimic the onset of failure, we take snap shots at
for the moderately elliptical cross-section and at
for the severely elliptical cross-section. These correspond to the first points of contact of the six different cross-sections of each set.
Figure 7 and
Figure 8 illustrate the evolution of the different failure modes for varying relative mucosal thickness ratios
and for different first points of contact
ϑ. Similar to
Figure 6, the upper right section displays the circumferential stress; the lower left section displays the three-dimensions of the failure mode. Again, red circumferential stress values at the inner mucosa correspond to high traction and blue values at the outer submucosa correspond to low traction. In elliptical cross-sections, the curvature is highest towards the long ends of the semi-major axis. This is where the initial instability occurs naturally. As growth progresses, the instability gradually propagates inward along the semi-major axis and forms additional folds until it reaches the semi-minor axis. As growth progresses further, the folds become more and more pronounced until first contact occurs. For regular ellipses, first contact always occurs at the long ends of the semi-major axis. Because of their inherent geometric heterogeneity, elliptical cross-sections tend to fold naturally in a double-symmetric folding pattern with naturally occurring outward folds along their semi-major axis. This implies that unperturbed simulations of elliptical cross-sections only generate even failure modes,
. As the relative mucosal thickness decreases, the number of folds increases from
at
and
to
at
and
, illustrated in the top and bottom rows of
Figure 7 and
Figure 8. In agreement with
Section 3.2, cross-sections with a relatively thin mucosal layer, towards the right column and the bottom row, require more growth to form the first contact. As the ellipticity increases,
the underlying failure mode at the first contact becomes more localized. Accordingly, cross-sections with larger ellipticity, from
Figure 6 through
Figure 8, require less growth to form the first contact. However, in agreement with
Section 3.1, since contact occurs only locally in cross-sections with larger ellipticity, the lumen area at the first contact remains significantly larger. This might indicate that circular cross-sections obstruct more drastically than elliptical cross-sections.
Figure 7.
Sensitivity of the failure mode with respect to relative mucosal thickness for moderately elliptical cross-section with . The instability occurs naturally at regions of highest curvature along the semi-major axis and propagates inward until it reaches the semi-minor axis. First contact occurs at regions of highest curvature and induces a moderately localized failure mode. This might indicate that airway obstruction is less drastic in moderately elliptical cross-sections than in circular cross-sections.
Figure 7.
Sensitivity of the failure mode with respect to relative mucosal thickness for moderately elliptical cross-section with . The instability occurs naturally at regions of highest curvature along the semi-major axis and propagates inward until it reaches the semi-minor axis. First contact occurs at regions of highest curvature and induces a moderately localized failure mode. This might indicate that airway obstruction is less drastic in moderately elliptical cross-sections than in circular cross-sections.
Figure 8.
Sensitivity of the failure mode with respect to relative mucosal thickness for severely elliptical cross-section with . The instability occurs naturally at regions of highest curvature along the semi-major axis and propagates inward until it reaches the semi-minor axis. First contact occurs at regions of highest curvature and induces a severely localized failure mode. This might indicate that airway obstruction is less drastic in severely elliptical cross-sections than in circular and moderately elliptical cross-sections.
Figure 8.
Sensitivity of the failure mode with respect to relative mucosal thickness for severely elliptical cross-section with . The instability occurs naturally at regions of highest curvature along the semi-major axis and propagates inward until it reaches the semi-minor axis. First contact occurs at regions of highest curvature and induces a severely localized failure mode. This might indicate that airway obstruction is less drastic in severely elliptical cross-sections than in circular and moderately elliptical cross-sections.
Figure 9 summarizes the above observations in two graphs to quantify the number of folds and the critical growth for varying ellipticity and varying relative mucosal thickness. With increasing relative mucosal thickness, the number of folds and the critical growth at the first contact decrease. With increasing ellipticity, the number of folds increases, while the critical growth at the first contact decreases. In view of the folded elliptical configurations displayed in
Figure 5,
Figure 7 and
Figure 8, this decrease of critical growth seems less alarming, since contact occurs only locally at the two outermost folds, while the lumen area remains relatively large. Overall, the two graphs confirm that smaller airways are at a higher risk of airflow obstruction than larger airways and that circular cross-sections obstruct more drastically than elliptical cross-sections.
Figure 9.
Number of folds and critical growth for varying ellipticities and varying relative mucosal thicknesses. With increasing relative mucosal thickness, the number of folds and the critical growth at the first contact decrease. With increasing ellipticity, the number of folds increases, while the critical growth decreases. This might indicate that smaller airways are at a higher risk of airflow obstruction than larger airways and that circular cross-sections obstruct more drastically than elliptical cross-sections.
Figure 9.
Number of folds and critical growth for varying ellipticities and varying relative mucosal thicknesses. With increasing relative mucosal thickness, the number of folds and the critical growth at the first contact decrease. With increasing ellipticity, the number of folds increases, while the critical growth decreases. This might indicate that smaller airways are at a higher risk of airflow obstruction than larger airways and that circular cross-sections obstruct more drastically than elliptical cross-sections.