## Abstract

The outward tethering forces exerted by the lung parenchyma on the airways embedded within it are potent modulators of the ability of the airway smooth muscle to shorten. Much of our understanding of these tethering forces is based on treating the parenchyma as an elastic continuum; yet, on a small enough scale, the lung parenchyma in two dimensions would seem to be more appropriately described as a discrete spring network. We therefore compared how the forces and displacements in the parenchyma surrounding a contracting airway are predicted to differ depending on whether the parenchyma is modeled as an elastic continuum or as a spring network. When the springs were arranged hexagonally to represent alveolar walls, the predicted parenchymal stresses and displacements propagated substantially farther away from the airway than when the springs were arranged in a triangular pattern or when the parenchyma was modeled as a continuum. Thus, to the extent that the parenchyma in vivo behaves as a hexagonal spring network, our results suggest that the range of interdependence forces due to airway contraction may have a greater influence than was previously thought.

- bronchoconstriction
- airway responsiveness
- computational model
- finite-element method

the lung parenchyma exerts outward tethering forces on the airways embedded within it. With lung inflation, these forces increase and, thus, progressively limit the ability of the airway smooth muscle to shorten, as demonstrated by the extreme sensitivity of airway responsiveness to transpulmonary pressure (6, 11, 14). This has led to the suggestion that a mechanical decoupling of the parenchyma from the airways might play an important role in the hyperresponsiveness of asthma (2). Local parenchymal distortion around a contracted airway might also influence the ability of another nearby airway to narrow.

Despite the clear physiological importance of airway-parenchymal interdependence, however, we have yet to understand the precise way in which it determines the hyperresponsive phenotype. In particular, our understanding of the magnitudes and topographic distributions of interdependence forces is based largely on theoretical results derived from a continuum mechanics analysis; yet, on a small enough scale, the lung is clearly not a continuum. For example, a very distal small airway may have only a handful of parenchymal attachments tethering its outer wall, a situation that would seem to be more accurately described by a discrete spring network, rather than a continuum. Indeed, two-dimensional spring network models have been used to model aspects of parenchymal function, such as the progression of emphysema (4, 25) and pulmonary fibrosis (4). Thus, while continuum theory has become widely accepted as applicable to airway-parenchymal interdependence, how it compares with models composed of networks of discrete springs is unclear. This is an important issue, given the key role played by small airway narrowing in the phenotype of airway hyperresponsiveness (2, 27).

To improve our understanding of this issue, we developed a numerical model of an elastic airway embedded in parenchyma. This model predicts the radial distributions of parenchymal distortions and forces that are produced when the airway contracts. In the present study, we examine how these distributions differ when the parenchyma is modeled as an elastic continuum vs. as a hexagonal spring network (HSN).

## METHODS

### Computational Modeling

#### HSN model.

We created a HSN model as follows. The parenchyma was represented as a square (dimensions *L* × *L*) planar network of identical hexagonally arranged springs. The hexagons were assumed to represent cross sections of individual alveoli. The positions of the network nodes (where the springs terminated) were fixed around the boundary of the network. An airway was created in the center of the sheet by designation of an approximately circular loop of springs as comprising the airway wall. The network nodes and edges within this loop were removed to create the airway lumen. We then assumed that each spring in the model had a smaller zero-load length than its current length and solved for static force equilibrium over the entire system to obtain the baseline configuration of the model. This produced an airway that was close to circular, with baseline (precontraction) diameter *D* (Fig. 1). In this way, we established a baseline configuration for the model that reflected the presence of prestress in all the alveolar walls and in the airway wall (although we found that the degree of prestress did not affect the results we present below).

The equilibrium configuration of the spring network was determined as follows. We began by assigning each spring a linear force-length (F-*x*) relationship (where *x* is change in length from zero-load state), such that its stiffness (*E*) was a constant given by
*k* and *b* are constants. In *Eq. 2*, we set the ratio *b/k* to be 10, in accordance with values estimated for alveolar wall material (10). The value of *k* for the springs in the airway wall was set arbitrarily to unity, and the remaining parameters were chosen to give the appropriate qualitative model behavior. Thus, *k* was 100 times greater for the airway wall than for the parenchyma, although we found that the results of our simulations were essentially unaffected by the relative stiffness of the airway wall and parenchyma.

To calculate what happens to the model following airway contraction, we followed the common practice, as described previously (15), of calculating the additional stresses of deformation caused by contraction under the assumption that the initial network configuration is unstressed (which gives the same normalized stresses and strains observed in the presence of prestress when the elastic elements are linear, because the stresses and strains in the structure follow the principle of superposition). A hexagonally arranged pin-jointed network of springs is not mechanically stable without prestress in the springs. In the present study, the equilibrium configuration of this structure was solved through a successive relaxation procedure using custom-purposed software based on standard finite-element methods (31) in which the structure was assumed to be rigid at each step in the process. That is, the edges were assumed to be Euler-Bernoulli beams interconnected via rigid joints where bending rigidity exists. The spring forces in each beam, according to their respective lengths, were then added to provide a resultant force at each link. These resultant forces were regarded as being externally applied to the frame structure, allowing a new network configuration to be derived on the basis of static force equilibrium. In computing this new geometry, the torsions at each link were ignored, so that the springs remained straight. Spring forces were recomputed based on the new edge lengths, and new resultant forces at links were obtained. This procedure was repeated until the maximum value of the resultant forces at all links throughout the network was smaller than the linear spring force at 0.1% extension.

Airway contraction was simulated by decreasing the zero-load length of the springs defining the airway wall. The zero-load lengths were decreased in small steps until a target reduction in airway area was achieved.

#### Triangular spring network model.

As a contrast to the hexagonal network, which is structurally unstable in the absence of prestress, we also modeled the parenchyma as a triangular network of springs having identical spring constants, as described previously (16). A triangular network of springs is stable in the absence of prestress. The network we used (Fig. 1*B*) had the same dimensions (*L* and *D*) as the HSN model. The springs in the triangular network were defined by *Eqs. 1* and *2* for the linear and nonlinear cases, respectively, and airway contraction was induced by application of inward radial forces to each of the airway wall nodes. The finite-element method was used to calculate the equilibrium configuration of the model, where each spring was represented as a linear bar element.

#### Continuum model.

Finally, to link our simulations to existing work on the modeling of airway-parenchymal interactions, we modeled the parenchyma as a linear continuum having the same dimensions (*L* and *D*) as the two network models described above. The model consisted of a triangular mesh (geometrically the same as Fig. 1*B*), with each triangular element being a linearly elastic constant strain element under the assumption of plane stress. Various values of parenchymal shear modulus have been reported in the literature, so we used a value of 5% of the bulk modulus based on measurements made in rat lungs (20). The exact value is not important for our purposes, however, because we are interested in the normalized displacement of the parenchyma around the airway for a given degree of airway narrowing, and an increase in shear modulus can be compensated for by an increase in the force of airway contraction to produce the same parenchymal displacement field. Airway contraction was achieved by application of an inward radial force on the nodes lying along the airway wall and calculation of the resulting force equilibrium configuration of the model using the finite-element method.

We define *r* to be the magnitude of the vector originating at the airway centroid and extending to a given point in the parenchyma after airway contraction. The normalized radial distance from this parenchymal point to the airway center is defined as *r*/*D*. After contraction, the parenchymal point moves a distance Δ*r* from its precontraction location. The normalized displacement of the parenchyma caused by airway contraction is then defined as the ratio Δ*r*/Δ*r*_{max}, where Δ*r*_{max} is the average displacement at the airway wall. On the basis of a theoretical analysis by Wilson (28), we know that the relationship between Δ*r*/Δ*r*_{max} = *X* and *r*/*D* = *Y* has the form
*A* and *B* are constants that depend on relaxed hole radius (*r*_{1}) and the ratio between *r*_{1} and the radius of the outer fixed boundary (*r*_{2}). *Equation 3* thus provides a means for testing how well continuum theory accounts for normalized radial displacement data as a function of *r*.

## RESULTS

### Parenchymal Force (Stress) Fields

Figure 2 shows the parenchymal forces generated by airway contraction for the HSN model (force in a random selection of springs) and for the continuum model (a single value of von Mises stress at each value of *r*). Force (F) is normalized to its maximum value (F_{max}) and is shown vs. *r*/*D* for the HSN and continuum models. Predictions for two different values of *L*/*D* are plotted for both models (the *L/D* values differed by a factor of 3, while F_{max} in the larger model was a factor of 0.128 and 1.02 of its value in the smaller model in the HSN model and the continuum model, respectively). Of particular note is the difference in the dependence of force (stress) on radial distance between the two models. Stress in the continuum falls off asymptotically toward zero and becomes negligible by about three or four diameters from the airway centroid. By contrast, the stress field for the HSN model changes rapidly within about three airway diameters (within this range, there are two populations of springs: one exhibits a decreasing stress with *r*/*D*, and the other an increasing stress). When *r*/*D* reaches a value of ∼3, however, stress in the HSN model reaches a plateau that extends to the parenchymal boundary. This plateau becomes more extended as *L*/*D* increases but decreases in magnitude, such that the total strain energy in the model (which is provided by the contraction of the airway) remains more or less unchanged.

In Fig. 3, we show that making the springs nonlinear in the HSN model, according to *Eq. 2*, has no important effect on the normalized displacement field or the normalized force field as functions of radial distance from the airway centroid.

### Parenchymal Displacement Fields

We contracted the airways in the three different parenchymal models (HSN, triangular spring network, and continuum) to achieve a 69% reduction in lumen area, as typically observed in experimental data from lung slices (1, 3, 12). Figure 4 shows normalized displacement vs. normalized radial distance for the HSN and continuum models. The HSN model predicted that displacement would fall off more slowly with distance from the airway than for the continuum model. The continuum model behaved virtually identically to the triangular spring network model (data not shown).

Figure 5 shows the fits provided by *Eq. 3* to the predictions of the HSN and continuum models. Points on and adjacent to the airway wall and the model boundary were ignored to avoid boundary effects on the fitted curves. Both gave good fits to the data, although the fit parameters were somewhat different (Table 1).

## DISCUSSION

Mead et al. (17) pointed out decades ago that the forces of mechanical interdependence between the intrapulmonary airways and the parenchyma surrounding them are of crucial importance to airway function. To understand these forces, they modeled the parenchyma as a two-dimensional hexagonal network of springs. Subsequently, the HSN model has been frequently invoked to explain important mechanical properties of the parenchyma, such as its stability (10, 24), and as a basis for modeling the progression of parenchymal diseases (4, 25). The aspect of airway-parenchymal interdependence that has excited the most consistent interest, however, relates to its putative role in limiting the ability of the airway smooth muscle to shorten. Dramatic evidence for this role is seen in the exquisite sensitivity of airway responsiveness to transpulmonary pressure. However, the extent to which the parenchyma actually affects airway contraction in healthy and diseased lung (18, 21) and whether abnormalities of airway-interdependence play a significant role in pathophysiology of asthma remain controversial (13, 27).

Because of the seminal work of Lai-Fook (15), most modeling studies of airway-parenchymal interdependence have considered the parenchyma to behave as a uniform and isotropic elastic continuum. This leads to a convenient expression for the way in which the outward tethering forces on an airway increase when the airway narrows concentrically (15), allowing the effects of transpulmonary pressure on airway responsiveness to be modeled in a manner that closely matches experimental observation (5). However, the interface between the parenchyma and the airway is not continuous but, rather, is mediated by a discrete number of alveolar septa. This raises questions about whether the parenchyma really does behave like an elastic continuum, especially when an airway is very small and may have only a few alveolar attachments around its circumference. Wilson (28) showed that a hexagonal network of springs can be approximated by a continuum mechanics model when parenchymal deformations are small and linear. However, when deformations are large around a contracting airway, the situation remains unclear.

Accordingly, we set out in the present study to compare the HSN and continuum models in terms of their predictions of stress and displacement fields radiating from a contracting airway. We found that the HSN and continuum models behave much as Wilson (28) predicted, with *Eq. 3* giving a good fit to the predictions of the radial parenchymal displacement field (Fig. 5). However, while the qualitative agreement between the two models is good in this sense, they differ rather substantially in their quantitative predictions (Fig. 4). Specifically, the HSN model predicts that the displacements in the parenchyma propagate considerably farther from the airway than does the continuum model. Furthermore, a triangular network of springs behaved similarly to the continuum model, so the differences in the model behaviors shown in Fig. 5 cannot simply be due to the discrete nature of the HSN model vs. the continuous nature of the continuum model. We speculate that these differences in behavior are more likely due to differences in structural stability. That is, a hexagonal network of springs can undergo shear deformation much more easily than can a triangular network because of the absence of diagonal struts in the HSN.

The above reasoning implies that the mechanical energy imposed by a contracting airway cannot be absorbed as readily through local network strains when the springs are arranged hexagonally, causing the strains to propagate farther from the airway than when the springs are arranged in a triangular pattern. Of course, this presupposes that the alveolar walls have negligible bending stiffness about their points of intersection, which seems reasonable if one assumes that the walls themselves are only capable of supporting linear tension. This supposition may be open to question, as some evidence of bending stiffness due to glycosaminoglycans has been reported (10). The long-range displacement propagation in a hexagonal network may be due, at least in part, to a mechanism known as isostaticity, in which the elastic equilibrium of a network can be obtained by considering only the static force equilibrium of the members impinging on each node, without the need to invoke the elastic properties of the individual members themselves (8). Isostaticity is likely not the only mechanism involved here, however, because force equilibrium at each node in our spring networks also involves changes in elastic energy due to changes in the lengths of the springs. Also, close to the airway wall, the parenchymal structure in the model (Fig. 1*A*) includes some pentagons (to accommodate the shape of the airway wall), so the isostaticity expected of a hexagonal network may have broken down to some extent in this region (8).

While the predictions of displacement field predicted by the HSN and continuum models differ only in a quantitative sense (Fig. 5), the predictions of the force (stress) field are different quantitatively and qualitatively (Fig. 2). The continuum model predicts a rapid and smooth monotonic decrease in stress that asymptotes toward zero, becoming negligible by about three or four diameters from the airway center. Furthermore, provided the parenchymal tissue is extensive enough to allow the stress to approach zero, its dependence on distance is independent of the size of the model. (The curves in Fig. 2 show the von Mises stresses, which we take as a key marker of overall tissue stress; the stress in a distorted continuum varies with direction, so it is impossible to capture every aspect of it in a single force-distance relationship.)

The HSN model also predicts that force near the airway changes rapidly with distance, but it is apparent that the springs near the airway wall break into two groups (we are able to readily show the force in each spring in this model). One group of springs is oriented approximately radially and is highly extended and, therefore, has high force, while the other group is oriented approximately circumferentially and has a much lower extension. By about three airway diameters from the wall, however, the spring forces converge toward a common value that remains at a fixed level virtually all the way to the fixed tissue boundary. As the boundary moves farther from the airway wall (i.e., as the ratio *L*/*D* increases), the plateau in force extends accordingly (Fig. 2). However, the magnitude of the plateau force decreases, so that the total strain energy in the model matches the energy put into contracting the airway. In fact, the stress normalization factor used for the larger HSN model in creating Fig. 2 is smaller than that used for the smaller model.

The HSN and continuum models thus make very different predictions as to how far elastic stress will propagate from a contracting airway. This has major implications for the role of airway-parenchymal interdependence in determining how the narrowing of one airway might influence the ability of another airway in its vicinity to contract. The continuum model predicts that airways will contract independently of each other, provided they are separated by more than about three precontraction diameters. The HSN model, on the other hand, predicts that the airway will always have some influence on each other's mechanical environment, albeit one that decreases progressively as they become more sparsely distributed throughout the parenchyma. It remains to be seen which model gives the better description of reality in the lung in vivo, but the matter has an important bearing on the way in which regional heterogeneities in bronchoconstriction develop. Recently, modeling work by Winkler and colleagues (26, 29, 30) has led to the proposition that the topographical distribution of ventilation defects in the constricted lung can be ascribed to a form of interdependence mediated by the constraint that a given volume of gas must be shared between a fixed number of airways. If the HSN model has any basis in reality, then our results suggest that airway interdependence may also be mediated to an important extent via the parenchyma in which the airways are embedded.

The HSN model predictions of stress propagation also have some interesting implications for the extent to which bronchoconstriction can increase overall parenchyma stiffness. Smith et al. (22) estimated that the airway tree contributes ≤20% to the recoil pressure of the lung and, therefore, calculated a limited influence of bronchoconstriction on overall lung stiffness. However, it has been shown in dogs (6) that when lung volume is maintained fixed, an intravenous injection of smooth muscle agonist can cause an 80–100% increase in static elastic recoil of the lung. A substantial stiffening of the parenchyma during airway contraction in rabbit has also been reported (19). These results thus support the notion that the stresses due to parenchymal distortion around contracting airways can be substantial in aggregate.

Our model results so far have been based on an assumption of linearity in the micromechanical behavior of the parenchyma. This makes for simplicity in the modeling process, but it is untenable as a description of reality. In fact, the one-dimensional stress-strain behavior of the parenchyma exhibits marked strain stiffening (7), so the individual springs in the HSN model ought to have appropriate nonlinear stress-stress relationships. Indeed, this is necessary to prevent the bulk modulus of the HSN model from decreasing as the model expands (23). We thus investigated how our above-described results are affected by making the HSN springs nonlinear according to *Eq. 2*. Figure 3 shows that nonlinearity has little effect on the shape of the normalized force and displacement fields. Of course, making the springs nonlinear can have a substantial effect on the overall magnitudes of the forces, but it is clear that our conclusions about propagation distances and the role of the parenchyma in mediating the energy of airway narrowing from one airway to another are not changed.

We have thus shown some significant, and perhaps somewhat surprising, differences between the predictions about airway-parenchymal interdependence that are made if the parenchyma is assumed to behave like a discrete network of hexagonally arranged springs as opposed to an elastic continuum. The predictions rest, of course, on a number of key assumptions that must be borne in mind. One of the most important is that we have modeled the parenchyma in only two dimensions. The choice of a hexagonal array of springs was motivated by the notion that this is a space-filling arrangement in two dimensions and that thin slices of lung parenchyma have a somewhat similar appearance under the microscope (27). The alveolar parenchyma, however, is essentially a three-dimensional foam, so it would have been more accurate to represent the parenchyma in two dimensions as a slice through a three-dimensional collection of packed dodecahedrons. Even that, however, neglects the fact that the alveoli are not all the same size and that the parenchyma contains numerous other nonairway structures, such as alveolar ducts and blood vessels.

Another issue for our model concerns the role of prestress. Although we arrived at the initial configurations of the network models by invoking prestress in the springs, we calculated the additional stresses due to airway contraction (Figs. 2 and 3) with the assumption that the initial conditions are stress-free. For linear springs, which obey the principle of superposition, the normalized force field due to airway contraction is independent of prestress. For nonlinear springs this is only approximately true. Nevertheless, to first order even for nonlinear springs, our results are independent of prestress and, therefore, have wide generality.

We could list numerous other approximations that might conceivably have an important influence on how well our model represents reality, such as not specifically including surface tension forces in the spring force-length relationships, having the parenchyma be purely elastic as opposed to viscoelastic, and assuming circular symmetry of the airway. Perhaps most important is our assumption that the parenchyma is mechanically uniform and isotropic. This is clearly not the case, as evidenced by a previous analysis of parenchymal strain fields surrounding contracting airways in lung slices (1, 9). Nevertheless, none of these factors is likely to change our overall conclusion that the structurally motivated HSN model and the conventional continuum model exhibit substantial differences in their mechanical behavior. These differences have important implications for the way in which a contracting airway influences the narrowing characteristics of other airways in the vicinity.

In conclusion, we have shown that when the parenchyma surrounding a contracting airway is modeled as a network of hexagonally arranged springs, the parenchymal stresses and displacements that are produced propagate substantially farther from the airway than if the parenchyma is assumed to be an elastic continuum. To the extent that the parenchyma in vivo behaves as a HSN, this suggests that the range of interdependence forces due to airway contraction has a greater influence than previously thought.

## GRANTS

This study was supported by National Heart, Lung, and Blood Institute Grant R01 HL-103405 and National Center for Research Resources Center of Biomedical Research Excellence Grant P20 RR-15557.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

B.M. performed the experiments; B.M. analyzed the data; B.M. and J.H.B. prepared the figures; B.M. and J.H.B. drafted the manuscript; B.M. and J.H.B. edited and revised the manuscript; J.H.B. are responsible for conception and design of the research; J.H.B. approved the final version of the manuscript.

## ACKNOWLEDGMENTS

The authors acknowledge helpful discussions with Michael J. Sanderson (University of Massachusetts Medical School).

- Copyright © 2012 the American Physiological Society