## Abstract

The forces of mechanical interdependence between the airways and the parenchyma in the lung are powerful modulators of airways responsiveness. Little is known, however, about the extent to which adjacent airways affect each other's ability to narrow due to distortional forces generated within the intervening parenchyma. We developed a two-dimensional computational model of two airways embedded in parenchyma. The parenchyma itself was modeled in three ways: *1*) as a network of hexagonally arranged springs, *2*) as a network of triangularly arranged springs, and *3*) as an elastic continuum. In all cases, we determined how the narrowing of one airway was affected when the other airway was relaxed vs. when it narrowed to the same extent as the first airway. For the continuum and triangular network models, interactions between airways were negligible unless the airways lay within about two relaxed diameters of each other, but even at this distance the interactions were small. By contrast, the hexagonal spring network model predicted that airway-airway interactions mediated by the parenchyma can be substantial for any degree of airway separation at intermediate values of airway contraction forces. Evidence to date suggests that the parenchyma may be better represented by the continuum model, which suggests that the parenchyma does not mediate significant interactions between narrowing airways.

- airway mechanics
- parenchymal mechanics
- finite element model
- spring network model
- mechanical heterogeneity

the airways and parenchymal tissues of the lung are intricately interwoven and thus tightly coupled to each other mechanically. This mechanical coupling, termed airway-parenchymal interdependence, has been extensively studied both experimentally and computationally (7, 17, 24, 25, 33) because the parenchymal attachments to the outside of an airway transmit transpulmonary pressure to the airway wall. This opposes the shortening of airway smooth muscle and so causes a strong inverse dependence of airway responsiveness on lung volume (6, 9). By the same token, when an airway narrows, it distorts the surrounding parenchyma for some distance out from the airway wall (23). Any other airway within this parenchymal distortion field will then have its own forces of airway-parenchymal interdependence increased accordingly.

The parenchyma can thus mediate between the mutual narrowing ability of two airways lying in sufficiently close juxtaposition, and thus has the potential to influence regional heterogeneity of airway narrowing, which is a key feature of obstructive lung disease (31). Whether these interactions actually occur to a significant degree in the lung, however, is not well understood. In fact, the only investigations of this issue to date have focused on the mutual interactions between large airways and their adjacent blood vessels (18, 19). Significant airway-airway interactions will only occur if the parenchymal distortion field generated around a contracting airway reaches as far as the nearest neighboring airway. In the present study, we consider two airways embedded in elastic parenchyma, and investigate how close the airways have to be in order for the intervening parenchyma to transmit a significant mechanical influence from one to the other. We made these calculations using three different two-dimensional computational models of the parenchyma, one consisting of a hexagonal network of springs representing the individual alveolar walls, the second consisting of a triangular network of springs, and the third consisting of an elastic continuum (21).

## METHODS

#### Overall modeling approach.

We consider two circular airways embedded in elastic parenchyma, modeled in two-dimensional cross section. The airways have identical relaxed diameter *d*_{0}, and their centers are placed at symmetrical positions with respect to a vertical line passing through the center of the tissue (Fig. 1). The square borders of the parenchymal tissue have a length of approximately 35*d*_{0}, and are fixed in place. In order to gauge the extent to which the parenchyma mediates between the mutually narrowing ability of the two airways, we simulated the narrowing of one of the airways (*Airway 1*) when the other airway (*Airway 2*) remained relaxed, and compared this to the narrowing of *Airway 1* when *Airway 2* also contracted. We then determined how the constricted area of *Airway 1* changed as a function of the distance between the two airways both when *Airway 2* was relaxed and when it was also constricted. Simulations were performed in Matlab (ver. 2010, MathWorks) on a desktop computer.

#### Modeling the parenchyma as a hexagonal spring network.

We first considered a two-dimensional slice of parenchyma consisting of a network of discrete elastic elements connected together in a repeating hexagonal pattern (Fig. 1), each elastic element representing an individual alveolar wall as seen in histological sections of lung parenchyma (10). The elastic elements were connected via frictionless pin joints. This representation has a long history of use as a model of mechanical interdependence and parenchymal stability (25, 28), and we have used it recently in detailed investigations of airway-parenchymal interdependence (21–23).

Each spring in the network follows a nonlinear force-extension relationship: (1)
where *y* is change in length relative to zero-load length, and *k* and *b* are spring constants. The precontraction model configuration corresponds to a prestressed state as illustrated in Fig. 1. This configuration was obtained by first creating a uniform spring network with each spring under tension (i.e., its zero-load length was less than its actual length) and then removing springs to create holes representing the two airways. Force equilibrium for the entire system was determined, and the resulting network configuration was used as the baseline for subsequent simulations of airway narrowing.

The airway wall was represented by a ring of springs bordering the hole representing the airway lumen. The wall of an airway is much thicker and stiffer than an alveolar wall (6), and the airway has some ability to resist changes in its lumen shape. To account for this, we incorporated a shape-maintaining mechanism by adding a radial force (*F*_{i}) acting at each spring junction around the wall in order to resist changes in airway shape according to (2)
where *i* denotes the *i*^{th} spring junction, *S*_{w} is a scale factor, *R*_{i} and *R*_{i0} are measures of lumen shape at the current and precontraction configurations of the model, respectively. *R*_{i} is the relative radial distance from the centroid of the airway lumen, (3)
where *r*_{i} is the radial distance of the *i*^{th} wall spring junction from the lumen centroid, and *r*_{m} is the average radial distance for all junctions. *F*_{i} defined in this way is sensitive to changes in airway lumen shape but not to changes in area.

As airway contraction occurs on a time scale of seconds to minutes (2, 11, 12, 16, 23), the deformation of the model due to airway contraction was regarded as a quasi-static process, and was calculated in a sequence of small time steps to produce a time course of change in airway area. Dynamic airway narrowing was achieved by having each wall spring follow the Hill hyperbolic relationship between force (*F*) and velocity (*v*) (6, 11), thus: (4)
where *F*_{0} is isometric force, *a* and *b* are model constants, *a* = *F*_{0}/4. For each time step in the simulation of contraction, a target change in the length of each spring comprising the airway wall was determined from the contraction velocity specified in *Eq. 4*. A tentative change in zero-load length was applied to the springs comprising airway wall and a new model configuration was obtained by solving for force equilibrium (described below). Based on the new configuration, the actual change in spring length was compared with the target value, a new value of zero-load length was tried, and the procedure was repeated until the target length was achieved for each spring.

To find the static elastic equilibrium configuration of the model at each time step, we used a successive relaxation approach based on the finite element method, following a strategy used to solve thermal stress in structures (36). To deal with the numerical difficulties presented by the fact that a two-dimensional pin-jointed hexagonal network is unstable without prestress (27, 28), we used the following three-step iterative procedure.

*Step 1*: Each spring is regarded as an Euler-Bernoulli beam with circular cross section with fixed connections between beams so that the whole network deforms as a statically determinate structure, based on a standard beam equation and stiffness equation for a plane frame element (36). The deformation of the network is determined given the incremental changes in the lengths of the airway walls (due to active smooth muscle contraction) and the shape restoring forces *F*_{i} applied to the airway spring junctions, as well as resultant nodal forces (which were zero for the first iteration).

*Step 2*: Update the model geometry by assuming all springs remain straight, and compute the resulting force in each spring due to the difference between its current length and the zero-load length. Compute the resultant force at each node in the network by vector summation of the forces in all connected springs.

*Step 3*: Using the resultant nodal forces computed in *Step 2*, repeat from *Step 1* until the nodal forces are negligible.

We used the following parameter values to solve the model, which consisted of 81 columns of hexagons each having an initial edge length of unity. For the purposes of determining the baseline geometry of the model, the zero-load length of the spring representing the edge of each hexagon was 0.9. The elastic constants for the parenchymal springs were *b* = 1 and *k* = 10 (in arbitrary units), while the corresponding constants for the airway springs were 100 times larger (6), mirroring a previous modeling study by Brook et al. (8). For the Euler-Bernoulli beam, the ratio of cross-sectional area to area moment of inertia was 100. An adaptive relaxation factor (*F*) was used to iteratively update nodal locations such that ** X** =

*X*_{0}+

*f*

**, where**

*dX***and**

*X*

*X*_{0}are the nodal position vectors at the current and previous iterations, respectively,

**is the predicted nodal displacement vector, and**

*dX**f*had an initial value of 0.5 and varied according to convergence speed. We allowed for a minimum lumen area at the end of contraction of 16% of baseline (before contraction), while the springs comprising the airway wall had a minimum allowed length of 40% of their length before contraction (3).

*S*

_{w}in

*Eq. 2*was set to 0.1. The values of

*b*and

*F*

_{0}were chosen by trial and error to obtain a desired degree of narrowing. The predicted contraction velocities thus were not calibrated to physiological values, and the forces generated by the smooth muscle in the airway wall are only meaningful relative to the stiffness of the springs in the parenchyma.

#### Triangular spring network model.

The motivation for using a hexagonal arrangement of springs in a network model is that this pattern represents, in a highly idealized manner, the space-filling appearance of adjacent alveoli in a thin slice of lung tissue. However, such an arrangement has an extremely low shear modulus because the hexagons offer essentially no resistance to shape change. In order to investigate a spring network model with a greater shear modulus than that of the hexagonal spring network model, we created a triangular network from the hexagonal network by connecting the six vertices of each hexagon by a spring to the center of the hexagon (Fig. 1). Each spring in the triangular network was identical and obeyed the same nonlinear force-length relationship as in the hexagonal spring network model. Airway contraction to 16% of the initial area was initiated by applying inward radial forces at the nodes around the borders of the two airways. Static force equilibrium was solved by applying the finite element method, assuming each spring to be a truss bar element.

#### Continuum mechanics model.

We also modeled the parenchyma as a linearly elastic continuum represented by a triangular mesh (Fig. 1). The mesh was denser toward the airway wall, with cells increasing in size by a factor of 1.2 away from the wall. The airways themselves were represented simply as circular holes in the continuum. The bulk modulus of the parenchyma was 1 and the shear modulus was 0.05. Airway contraction was achieved in a single step by applying equal radial forces inward at each node around the airway lumen. The force was chosen to achieve a contracted lumen area ≥16% of baseline. This degree of contraction was chosen to match that of the spring network models above for ease of comparison. Accordingly, the magnitude of the contraction force was thus not meant to be physiological; it has meaning only relative to the elastic moduli of the parenchyma. Similarly, the contraction forces used in each of the three models are not comparable; each model serves as its own control in terms of revealing how a given degree of luminal narrowing affects the distortion of the surrounding parenchyma and the radial propagation of parenchymal forces. The deformation of the model under airway contraction was determined using the finite element method (36), assuming a linear constant strain element at each triangle.

## RESULTS

In the hexagonal spring network model, the dynamics of airway narrowing were controlled by the values of the parameters *F*_{0} and *b*, which control contraction force and velocity, respectively. When both airways contracted with the same value of *F*_{0} and different values of *b* (0.4 for *Airway 2* and 0.2 for *Airway 1*), their contraction time histories were different, but they eventually achieved comparable lumen areas (Fig. 2). This demonstrates that the final degree of airway narrowing in this model is predominantly determined by smooth muscle force rather than shortening velocity.

The final degree of airway narrowing in the hexagonal spring network model was influenced by whether only one or both of the airways contracted, but only if the force of contraction was not so large as to completely overpower any elastic forces generated within the parenchyma. For example, with the contraction force produced by setting *F*_{0} = 0.005, the decrease in lumen area of *Airway 1* was about twice as great when *Airway 2* remained relaxed compared with when it also narrowed (Fig. 3*A*). By contrast, when the contraction force was large (*F*_{0} = 0.01), the activation state of *Airway 2* had no noticeable effect on the narrowing of *Airway 1* until two airways were within about three precontraction diameters of each other (Fig. 3*B*). Likewise, we also found that if the airway contraction forces were very small, the parenchymal distortion forces generated between them were also too small to have an influence. It is thus only at intermediate levels of contraction force that we find the contraction of one airway influences the narrowing of another nearby airway in this model. Accordingly, because the forces of airway interdependence increase with transpulmonary pressure, it is possible that the ability of one airway to influence another could depend markedly on lung volume.

We also note that there are two distinct regions in the plots of airway narrowing vs. separation in Fig. 3. When the airways were relatively far apart their degrees of airway narrowing did not change with separation distance. This is because of the stable plateau in parenchymal distortion force that arises further than about three airway diameters when an airway contracts within a hexagonal spring network, as we have shown in previous publications (21–23). We refer to this region as the far field. Conversely, when the airways were close together the degree of narrowing depended on separation distance. We refer to this region of the plot as the near field, and it is caused by the fact that as the airways approach each other, their immediate environments become increasingly dominated by the hole in the parenchyma that constitutes the lumen of the other airway. We found, however, that the relative influence of one contracting airway on the narrowing ability of another airway was not significantly affected by whether the two airways were within each other's near field or far field. Fig. 4 shows the relationship between the narrowing of *Airway 1* without vs. with the narrowing of *Airway 2* as a function of contraction force. The narrowing ratio peaks at a value somewhat over 2 when *F*_{0} is around 0.005 for a separation distance taken both from the near field and from the far field.

The above observations were qualitatively similar at other ratios of parenchymal network size to airway diameter.

For the continuum model, we found that the degree of narrowing of *Airway 1* is independent of the contraction of *Airway 2* for all values of contraction force except when the airway separation distance is less than two precontraction diameters, at which point the contraction of *Airway 2* created a very slight impediment to the narrowing of *Airway 1* (Fig. 5). Consequently, the relationships of *Airway 1* narrowing ratio vs. contraction force were virtually flat for both the near field and the far field (Fig. 6).

We also found that contracted area exhibited a positive dependence on separation distance when the two airways were close together in the continuum model (Fig. 5), in contradistinction to the behavior exhibited by the spring network model (Fig. 3). This opposite behavior of the two models is simply an artifact of the way their respective airways were represented. In the spring network model, the airway wall was stiffer than the parenchyma and had some degree of resistance to changing shape, so as the two airways approached each other they presented increasing impediments to each other's ability to narrow. By contrast, the airways in the continuum model were represented simply as holes in the parenchyma, so the proximity of one airway to another presented a reduction in impediment to narrowing.

For the triangular spring network model, the contracted area vs. separation distance relation (Fig. 7) closely follows that of the continuum model, suggesting minimal influence between the two airways when contracting, except when the airway separation distance is less than two precontraction diameters.

The above results apply to the case when the two airways are of identical size. Simulations were also performed for all three parenchymal models when *Airway 2* was twice as large as *Airway 1*. The results were very similar to the case of uniform airways, i.e., little interairway interaction for continuum and triangular spring network models, and conditional interaction for the hexagonal spring network model. Thus all conclusions in this study apply to airways of both uniform and nonuniform sizes.

## DISCUSSION

Regional heterogeneity of lung function is a general feature of normal as well as diseased lung (15). Furthermore, heterogeneity per se, independent of the mean level of airway smooth muscle activation, can have a major negative impact on overall lung function (5, 15, 31). The genesis of regional heterogeneity of airway narrowing in the lung is usually assigned to anatomic variations in effectors such as smooth muscle mass and airway geometry that, when modeled in combination, can give rise to wide distributions of airway luminal narrowing throughout the lung (14, 20, 29, 30, 32). Recently, heterogeneity of airway narrowing has also been linked to emergent phenomena that potentially arise from the interactions taking place between airways as they endeavor to share gas volume under the constraints imposed by their competing contraction dynamics (31, 34, 35). Such interactions, however, are mediated via the airway lumen as a result of the involved airways drawing their respective flows from a common flow pool. In the present study, we explore a different kind of airway-airway interaction, one that is mediated via the intervening parenchyma.

The results of our study indicate that if the parenchyma of the air-filled lung in vivo behaves like an elastic continuum, then contracting airways would have to be very close together (within about two airway diameters) to have even a slight influence on each other's ability to narrow (Figs. 5 and 6). Such a degree of closeness between airways seems unlikely based on the following order-of-magnitude argument. The airways and alveoli occupy, to a rough approximation, 2,000 ml vs. 100 ml, respectively. If one assumes that the alveoli and airways are scattered randomly throughout the volume occupied by the lung, then one would expect any subset of the lung volume to be occupied by airway lumen vs. alveolar lumen in the ratio 1/20 (assuming that the volume taken up by actual tissue can be ignored). A random cross section through this volume would thus consist of a patchwork of airway and alveolar lumens occupying areas in the ratio 1/20 raised to the two-thirds power, which is the ratio 1/7.4. Applying the same dimensional reasoning, a random line through such a cross section would traverse the lumens of alveoli and airways in the ratio 1/20 raised to the one-third power, which is 1/2.7. In other words, if we assume that the average diameter of the airways intersected by the line is 1 length unit, then on average these airways will be separated by 2.7 length units of alveolar parenchyma, which approaches the far-field regions of the airway interaction curves shown in Figs. 5 and 7 for the continuum and triangular spring network models, respectively.

On the other hand, if the parenchyma of the air-filled lung in vivo behaves like a two-dimensional hexagonal spring network, then the airways would be expected to significantly influence each other's ability to narrow at any separation distance provided that the forces of smooth muscle contraction do not overwhelm the distortional forces generated in the parenchyma, as illustrated in Figs. 3 and 4. Thus the potential for one contracting airway to influence the narrowing of another depends on the force exerted by the airway smooth muscle; too little and the parenchymal forces are inconsequential, too large and parenchymal forces are negligible in comparison. When such interactions do occur, however, we predict that they serve to reduce the regional heterogeneity of narrowing throughout the lung because parenchymal distortion causes mutual inhibition of narrowing among interacting airways, forcing them toward a common end point (Figs. 2–4). However, we do not yet know whether this actually applies to a real lung because the applicability of the hexagonal spring network model has not been established.

Continuum mechanics analysis (1, 8, 21, 33) gives rise to the conventional view that the parenchymal deformation field induced around a contracting airway reaches out to only about three diameters from the airway wall. This has been supported by mechanical tests on explanted lung slices (2, 8, 12, 16, 23), but these slices were filled with agarose to stabilize the otherwise highly collapsible alveolar structures, so they may have behaved more like a continuum than if they had been filled only with air (8, 21). Furthermore, we (21) recently found that if the parenchyma is modeled in two dimensions as a hexagonal spring network, the distortion forces induced by a contracting airway can propagate throughout the entire parenchymal domain.

Whether the parenchyma of the air-filled lung in vivo behaves more like a continuum or a hexagonal spring network remains an open question. The source of the difference in behavior between the two models seems to arise from their differences in shear modulus relative to bulk modulus (21). Therefore, to gain further insight we considered an intermediate case in which the very small shear modulus of the hexagonal network was increased to be similar to that of the continuum model by adding cross springs to each hexagon to create a triangular pattern as shown in Fig. 1. The resulting behavior (shown in Fig. 7) was similar to that of the continuum model (shown in Fig. 5), suggesting that if alveolar walls behave as elastic sheets that resist shape change to any significant degree then the parenchyma itself might be more likely to behave like the continuum model. Indeed, there is reason to suspect that this is the case because the surface forces in the air-liquid interface (13, 28) and the stabilizing effects of proteoglycans (10) could both act to stabilize alveolar wall shape. Nevertheless, this still does not prove the matter definitively one way or the other.

We also made other simplifying assumptions that could have affected our results. For example, we assumed the parenchyma to be homogeneous and isotropic, although we have recently demonstrated that parenchymal heterogeneities are unlikely to have a major effect except in rather special cases (22). Similarly, we assumed regional ventilation throughout the lung to be homogeneous, whereas in reality it can be quite variable resulting in different degrees of small airway inflation, so our results pertain only to an average airway. We also assumed that while the parenchyma affects the mechanical behavior of the airways, there is no influence in reverse. In fact, it has been shown that the airway narrowing can increase parenchymal stiffness (26), presumably because of the distortion caused in what is a nonlinearly elastic material. Finally, we assumed the stiffness of the alveolar wall to be 100 times that of the springs in the parenchyma. However, experimenting with stiffness ratios of 1 to 200 gave similar results, indicating that this ratio is unimportant for our overall conclusions.

These assumptions aside, however, our results indicate that if the parenchyma behaves like a mechanical continuum then significant airway-airway interactions mediated by the parenchyma during bronchoconstriction are highly unlikely, while such interactions could be significant if there is any validity to the hexagonal spring network representation of the parenchyma. In the former case, we would expect airway-airway interactions to have no influence on the regional heterogeneity of airway narrowing, while in the latter case we would expect heterogeneity to be reduced with increasing lung volume as the parenchymal mediation forces increase. There is some evidence that heterogeneity of narrowing increases with decreases in lung volume (4), but whether this is due to the mechanisms considered remains to be determined.

In summary, our results indicate that while the importance of airway-airway interactions mediated by the parenchyma during airway narrowing remains an open question, the weight of evidence at this point would suggest these effects to be very minor. This is good news for efforts to computationally model the mechanical behavior of the lung in an anatomically and physiologically faithful fashion because it obviates the need to take these interactions into account, something that so far has not been done, and which would be far from trivial.

## GRANTS

This work was supported by National Heart, Lung, and Blood Institute Grant R01 HL-103405 and National Institute of General Medical Sciences Grant P30 GM-103532.

## DISCLOSURES

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

## AUTHOR CONTRIBUTIONS

Author contributions: B.M. and J.H.B. conception and design of research; B.M. performed experiments; B.M. analyzed data; B.M. and J.H.B. interpreted results of experiments; B.M. and J.H.B. prepared figures; B.M. and J.H.B. drafted manuscript; B.M. and J.H.B. edited and revised manuscript; B.M. and J.H.B. approved final version of manuscript.

- Copyright © 2014 the American Physiological Society