Journal of Applied Physiology
Vol. 82, No. 1, pp. 32-41, January 1997
GAS EXCHANGE, MECHANICS, AND AIRWAYS

A distributed nonlinear model of lung tissue elasticity

Geoffrey N. Maksym and Jason H. T. Bates

Meakins-Christie Laboratories and Department of Biomedical Engineering, McGill University, Montreal, Quebec, Canada H2X 2P2

ABSTRACT

INTRODUCTION

METHODS

RESULTS

DISCUSSION

ACKNOWLEDGEMENTS

FOOTNOTES

REFERENCES

ABSTRACT

Maksym, Geoffrey N., and Jason H. T. Bates. A distributed nonlinear model of lung tissue elasticity. J. Appl. Physiol. 82(1): 32-41, 1997.- We present a theory relating the static stress-strain properties of lung tissue strips to the stress-bearing constituents, collagen and elastin. The fiber pair is modeled as a Hookean spring (elastin) in parallel with a nonlinear string element (collagen), which extends to a maximum stop length. Based on a series of fiber pairs, we develop both analytical and numerical models with distributed constituent properties that account for nonlinear tissue elasticity. The models were fit to measured stretched stress-strain curves of five uniaxially stretched tissue strips, each from a different dog lung. We found that the distributions of stop length and spring stiffness follow inverse power laws, and we hypothesize that this results from the complex fractal-like structure of the constituent fiber matrices in lung tissue. We applied the models to representative pressure-volume (PV) curves from patients with normal, emphysematous, and fibrotic lungs. The PV curves were fit to the equation V = A  Bexp(KP), where V is volume, P is transpulmonary pressure, and A, B, and K are constants. Our models lead to a possible mechanistic explanation of the shape factor K in terms of the structural organization of collagen and elastin fibers.

lung tissue mechanics; collagen; elastin; stress strain; power law; pressure-volume curves

INTRODUCTION

THE ELASTIC RECOIL of the lung at normal breathing frequencies is dominated by the nonlinear stress-strain characteristics of the lung tissue (35). It is commonly agreed that the main stress-bearing constituents of lung tissue are collagen and elastin fibers (30). These fibers differ significantly in their mechanical properties, with collagen being three to four orders of magnitude stiffer than elastin (40). Furthermore, elastin is stretchable to >250% of its original length before rupturing, whereas collagen can stretch merely 1-3%. The origin of the curvilinear stress-strain behavior is generally thought to be one of collagen fiber recruitment. That is, elastin is responsible for load bearing at low strains when much of the collagen is "wavy" and, therefore, not contributing to the tension. As strain increases, the collagen fibers become straight and so progressively take up more load, thereby stiffening the tissue. The amounts of collagen and elastin and their spatial arrangements in the lung are altered in diseases such as pulmonary emphysema and fibrosis (39, 43), and it is thought that these alterations are largely responsible for the abnormal pressure-volume (PV) curves associated with these diseases (27). However, it is currently unclear how to precisely translate alterations in lung tissue constituents into changes in tissue mechanical properties. We have developed a model of lung tissue that exploits the collagen fiber recruitment concept by representing the collagen and elastin fibers as a series of spring-string pairs. We explore how distributions of spring stiffnesses and string lengths can account for measured stress-strain curves of lung parenchyma. We also explore how our model might be used to explain the curvilinear quasistatic PV curve of whole lung.

METHODS

Analytic Model Development

Distributed stop length model. A preliminary development of this model applied to lung tissue strips was previously presented by the authors (17).

We consider the model of Fig. 1 with identical spring constants for each spring-string pair, whereas the stop lengths are distributed. At rest, there is no load on the model, and its resting length is L_r. When the model is extended to length L > L_r, all the units with stop lengths less some value l_o will have become stopped. The length of the model is then the sum of three terms

(1)

where L_s is the combined extensions of all the stopped units and L_u is the length contribution of the remaining units. The stop lengths (l) are distributed according to a density distribution N(l). Equation 1 becomes

(2)

As the model is lengthened, l_o increases and more units become stopped. Differentiating, we have

(3)

which gives the model compliance

(4)

where dl_o/dF is the compliance of the string-spring unit that has just become stopped and F is tension. We have defined this compliance to be the same for all units; thus

(5)

where k is the elastic constant. Differentiating Eq. 4 with respect to l_o gives

(6)

Recognizing that this is true for every value of l_o and substituting in for dl_o/dF from Eq. 5 we have

(7)

Equation 7 allows us to calculate N(l) to within the arbitrary elastic constant k from measurements of tissue force and length.

We note that Eq. 5 substituted into Eq. 4 and evaluated for l_o = 0 gives the initial slope of the length-tension curve, thus

(8)

where X is the total area under the stop length density distribution. We have no guarantee that X will be finite, so we merely ensure that the upper limit of the integral in Eq. 8 is at least as great as the maximum extension of the tissue generating the data to which Eq. 7 is applied. In this paper, we use an upper limit corresponding to a strain of 2 and a value for X of 10,000, thus defining a value for k.

Distributed stiffness model. An alternate model based on Fig. 1 can be developed in which the stiffnesses are distributed but the stop lengths are identical for all spring-string units. Let M(k) represent the stiffness density distribution. Equation 1 now becomes

(9)

where k_o is the stiffness of that unit which has just become stopped. Differentiating with respect to the force, we get the compliance of the model

(10)

Differentiating again with respect to k_o gives

(11)

Because all units of the model experience the same force, equal to the force of the unit that has just reached its stop length, we have

(12)

Consequently

(13)

Substituting Eq. 13 into Eq. 11 and recognizing that the result holds for all values of k_o gives

(14)

Equation 14 allows us to calculate M(k) to within the arbitrary stop length constant l from measurements of tissue force and length.

We note that the maximum extension (L_{s max}) is given by

(15)

where Y is the total area under the stiffness distribution. Similar to the case for the stop length distribution, we have no guarantee that Y will be finite, so we fix the upper limit of the integral in Eq. 15 such that the maximum strain is 2 when Y has a value of 10,000, thus determining l.

Curve Fitting to Experimental Data

(16)

To calculate the stop length distributions (Eq. 7) and stiffness distributions (Eq. 14) from stress-strain tissue data, the following relationship was used to convert from stress-strain to length-tension

(17)

where L_r and A_r are given for each tissue strip.

Computer Model Simulations

Table 1. Parameters from the tissue strips and the curve fits of the tissue strip stress-strain curves to Eq. 10, together with means and standard deviations Tissue Strip Ar, cm2 Lr, cm G, kPa H RMS residual, kPa 1 0.229 2.8 0.254 0.300 0.108 2 0.204 3.0 0.199 0.326 0.028 3 0.240 2.6 0.137 0.249 0.163 4 0.207 3.1 0.164 0.348 0.086 5 0.272 2.6 0.390 0.297 0.036 Mean ± SD 0.230 ± 0.025  2.8 ± 0.2  0.229 ± 0.100  0.304 ± 0.037  0.088 ± 0.055 RMS residual, root mean squared error from the fits. Ar, resting area; Lr, resting length; G and H, fitted constant.

RESULTS

Application to Lung Tissue Stress-Strain Curves

Applying the distributed stop length model via Eqs. 7 and 17 to the stress-strain Eq. 16, the stop length distribution is

(18)

where k is computed from Eq. 8 as

(19)

The resulting N(l) for each tissue strip are shown in Fig. 3, where the range of l over which N(l) is plotted corresponds to the measured forces from the relation

(20)

Similarly, applying the stiffness distribution model via Eqs. 12, 14, and 17 to Eq. 16 gives

(21)

where l is calculated from Eq. 15 as

(22)

where _f is the maximum strain of the model. In this case, k is our dependent variable and is chosen to correspond to the range of available data from the relationship

(23)

The resulting stiffness distributions for each of the tissue strips are shown in Fig. 4.

Computer Simulations

(24)

where F_f is the final force defined by _f = 2. In the distributed stiffness model, stiffnesses were distributed according to Eq. 21, with l from Eq. 22. The range of k was limited between 0 and

(25)

where again F_f is defined by _f = 2. The distributions for the combination model were created according to Eqs. 18 and 19, respectively, with parameters altered so as to more closely match the mean tissue stress-strain curve. Specifically, the value of k in Eq. 18 and in Eq. 24 decreased by a factor of 2 over the distributed l model, and the value of l of Eqs. 21 and 25 was increased by a factor of 2 over the distributed k model.

DISCUSSION

Tissue Models

The idea that elastin determines the stress-strain behavior of lung parenchyma at low strains, whereas collagen takes over as strain increases, is well known, first having been introduced by Setnikar (31) and extended later by Mead (20). More recently, there has been a few studies that explore this idea in a quantitative fashion. For example, Fung (10) used a linear expression for the constitutive equation for elastin together with an exponential function for collagen. The two fiber types were assumed to be distributed probabilistically according to observed distributions of fiber diameter and curvature. This gave a tissue strain energy function from which the incremental material constants were calculated according to the theory of continuum mechanics. Suwa et al. (37) ascribed the origin of the curvilinear stress-strain relationship of lung parenchymal strips solely to the reorientation of the constituent fibers, where a fiber does not support load until oriented in the direction of the applied strain. Models based on the gradual recruitment of kinked or "wavy" fibers were proposed by Soong and Huang (34) for human alveolar wall and by Decraemer et al. (6) for various biological soft tissues. However, these models consist of only one fiber type, elastic collagen, and assumed a specific functional form for the distribution of fiber lengths rather than estimating the distribution from the available data. Brown et al. (4) studied ligamentum propatagiale from birds and discussed a model similar to ours in which the collagen fibers had distributed initial waviness. However, they did not develop it mathematically and concluded that it may be an important factor in distensible tissues but probably was not the predominant mechanism in ligamentum propatagiale. Recently, Denny and Schroter (7) constructed a model for the alveolar duct including collagen and elastin fibers of three different diameters for different alveolar locations, based on the densities of collagen and elastin fibers found in dog lungs (25). They assumed a particular stress-strain relationship for the collagen fiber bundles, which they assumed arose from gradual recruitment of tensed collagen fibers within a bundle. Our approach, however, has a significant advantage over all the above studies in that it allows us to determine functional forms for the distributions of elastin or collagen fiber properties from experimental stress-strain data, rather than having to assume these forms a priori.

The distributions of elastin fiber diameters were measured in human pulmonary alveolar walls of human lungs by Sobin et al. (33) and in alveolar mouths and ducts by Matsuda et al. (19). The distribution of elastin diameters would be expected to relate to our distribution of elastin stiffnesses, given that the stiffness of a fiber or fiber bundle is proportional to its cross-sectional area for a fixed Young's modulus. Thus a plot of the stiffness distribution functions vs. the square root of stiffness should resemble a fiber diameter distribution function. This is borne out by the mean diameter distribution from our five tissue strips shown in Fig. 6, which demonstrates a striking similarity to the data of Sobin et al. (33) (Fig. 7). Distributions of collagen fiber curvatures in the lung (19, 33) are also similar in form to our stop length distributions, although it is not entirely clear how these two quantities relate.

These morphological observations lend support to the model, but this does not imply direct correspondence between the strings and springs of our model and the individual components of actual tissue. The arrangement of the fibers in lung tissue is not only serial but is also parallel and interwoven, forming complex interconnections (22). In fact, our stop length distribution may span multiple length scales. That is, the shorter collagen fibrils may correspond to the smallest stop lengths, but, as many fibers become taut, a stop length may be formed from a cluster of collagen fibers (connected in parallel and in series), all reorienting and straightening against the load. Thus the distribution of stop lengths in our model may represent how the collagen matrix is arranged in a fractal sense.

In our model, the properties of collagen and elastin are considered independently but their interaction is defined by their pairwise coupling and serial connection. Mercer and Crapo (22) showed that elastin and collagen fibers are interwoven especially in the alveolar entrance ring region, implying mechanical coupling of some kind. In fact, Mijailovich et al. (23, 24) considered this interaction to be responsible for the observed dynamic viscoelastic properties of lung parenchymal strips in a model based on sliding filaments. It may thus be that interaction along the lengths of the fibers is responsible for dynamic tissue properties, whereas interconnections between the two fiber matrices with fiber properties distributed in a fractal-like manner are responsible for the static mechanics.

We note that for the large majority of the data (>95%) the distributions of Figs. 3 and 4 are nearly power-law in form, manifest in the nearly linear tails, where N(l) is proportional to 1/l² and M(k) is proportional to 1/k. There are other examples of inverse power laws exhibited by the lung. For example, the stress response of lung tissue strips to a step change in strain exhibits power-law dependence with time. Also, when strips are mechanically cycled over a range of frequencies, the impedance amplitude is an inverse power-law function of the cycling frequency (3). The lung inflation process has also been found to follow power laws both as the magnitude of decreases in airway resistance measured from the terminal alveolar space and as the times between airway opening events (36). Complex systems with many varied elements interacting over a wide range of length scales tend to exhibit some measure or multiple measures in an inverse power-law sense of the form 1/f , with 0 <  < 2, and f being usually a distributed frequency, amplitude, or spatial measure. In particular, the inverse power-law function often arises as the result of a cascade of smaller events that occur in a chain reaction (2, 41). We hypothesize, therefore, that the bulk stress-strain properties of lung tissue are due in large part to some kind of fractal-like structure arising from biological self-organization. This is responsible for the 1/f character found for our collagen and elastin distributions, which is a manifestation of general complexity.

In each of our models, the parameter of the element in parallel with the distributed parameter is identical for the entire model. That is, for the distributed stop length model the stiffnesses are fixed, whereas in the distributed stiffness model the stop lengths are fixed. This is clearly not realistic considering natural physiological variability in real tissue. We thus wondered how robust our models are to variability in the fixed parameters of each model. We therefore generated an l-distributed model as before (Eq. 18) but with the stiffnesses chosen randomly from uniform distributions having mean k chosen as in Eq. 19 and with ranges from 0 to 100% of the mean. We found that the stress-strain curve changed only slightly, showing visible weakening from the mean tissue curve for ranges >50% of the mean (Fig. 8, long-dashed lines). Similarly, we generated a stiffness distributed model as before according to Eq. 21 but with uniformly distributed l of the same mean as the fixed l model and with increasing ranges from 0 to 100% of the mean. Again, this caused little change in the stress-strain curve, with some stiffening observed at ranges >50% of the mean (Fig. 8, dotted lines).

We see therefore that the model is robust to the introduction of variability in the fixed parameters. Nevertheless, it is obvious that both fiber types will be distributed in real tissue, as indicated by the data of Matsuda et al. (19) and Sobin et al. (33). We presumed that the forms of the collagen and elastin distributions are those given in Eqs. 18 and 21. However, each of these distributions accounts for the entire stress-strain characteristics independently. Thus, to invoke both simultaneously, we softened their individual contributions by doubling all stop lengths and halving all stiffnesses. This produces a stress-strain curve that matches the mean experimental curve fairly well, dissociating only at high strains (Fig. 5, dotted-dashed line).

Lung Tissue Diseases

We generate static elastic PV curves by having both stop volumes and unit elastances distributed according to Eqs. 7 and 14 with F, L, k, and l substituted for P, V, E and , respectively, where V is lung volume and P is transpulmonary pressure. The equation most commonly used for the static PV curve is that first introduced by Salazar and Knowles (28)

(26)

where A, B, and K are constants. If we solve for the model distributions by double differentiation of Eq. 26 in the PV equivalent of Eqs. 7 and 14, we have

(27)

where

(28)

and

(29)

where

(30)

and we have used the fact that

(31)

in analogy with Eq. 12.

Equations 27 and 28 represent the stop volume distribution model, whereas Eqs. 29 and 30 represent a distribution of elastances. It is likely that the collagen fiber matrix of lung tissue and the elastic constitutive elements have distributed stop volume and stiffness properties, respectively. We therefore recognize that the two models represent the extremes of possibility, and so consider these two extremes to quantitatively assess how changes in PV curves may reflect pathological changes in lung tissue.

We apply our model Eqs. 27-30 by using the values for the parameters of the Salazar and Knowles (28) equation (Eq. 26) obtained from the PV curves of representative emphysematous and fibrotic lungs (Table 2). The changes of the model parameters expressed as a percentage of normal are shown in Table 3. For emphysema,  and  are elevated, whereas E and  are decreased below normal. The stop volume distribution for an emphysematous lung is thus shifted toward large stop volumes compared with normal (Fig. 9). This is consistent with the nature of the disease, which consists of destruction of tissue with associated remodeling of both elastin and collagen networks as well as an increase in air-space volumes (12, 43). By contrast, in fibrosis, the changes in the tissue involve widespread inflammation, an increase in the concentration of interstitial fiber constituents, and air spaces are reduced or even filled in with inflammatory products (8). Correspondingly for fibrosis, parameters  and  are greatly decreased below normal while E and  are elevated, which results in reduction in vital capacity and a stiffening of the tissue. Note that  and  have units of volume (and are in fact proportional to B), and thus control the vital capacity of the model, whereas E and  have units of elastance and thus, with  and , affect the shape of the PV curve.

Table 2. Parameter values from the fit of Eq. 26 to the deflation limb of pressure-volume curves of a normal, emphysematous, and fibrotic lung in vivo A, liters B, liters K, cmH2O1 Normal 5.65 4.79 0.143 Emphysema 7.66 6.98 0.325 Fibrosis 2.64 0.93 0.089

Table 3. Percent changes in model parameters from normal values to those for emphysema and fibrosis for the distributed stop volume and distributed elastance models Parameter Distributed Stop Volumes, %  Distributed Elastances, %  E       Emphysema 30 146 30 146 Fibrosis 827 19 827 19 E, elastance.

A possible mechanism for the dominance of changes in elastance-related properties over changes in volume-related properties via our model is shown with the aid of Fig. 10. Here we consider the collagen and elastin fibers to be composed of the parallel contribution of multiple fibrils. The stiffness of a single unit is equal to the sum of the stiffnesses of all the individual elastin fibrils in that unit, whereas the stop length of the unit is the minimum collagen fibril length. A decrease in the number of elastin fibrils would manifest itself as a proportional decrease in the stiffness of the unit, whereas an increase in stop length would occur only if the shortest fibril were removed. The situation is similar for the addition of collagen or elastin fibrils: the elastance increases in proportion to the number of added elastin fibrils, whereas the stop length shortens only if a collagen fibril shorter than the shortest is added. Thus the elastance-associated parameters are more likely to be affected than volume-associated parameters under conditions that randomly affect both collagen and elastin, through tissue destruction (emphysema) or tissue building (fibrosis). Possible support for the model of Fig. 10 was demonstrated by Mercer and Crapo (22), who observed parallel collagen fibrils meandering in a wavelike fashion. They measured the maximum possible degree of fiber extension to be as much as 16.1 ± 3.2% (SE, n = 6) near the alveolar septal edge in human lungs. This is comparable in magnitude to the changes expected to occur between functional residual capacity and total lung capacity (21).

A significant result of the application of the model to the PV curves of diseased lungs is that it may provide a mechanistic description for the origin of the shape parameter K (from Eq. 26) in terms of the structural constituents of lung tissue and their organization. The parameter K has been found to be increased in those patients with emphysema (12, 26) and to be decreased in fibrosing alveolitis (38), as shown in Table 2. Within the context of our model, these changes in K are interpreted in terms of alterations to the collagen and elastin networks (Table 3). In fact, K is inversely proportional to the stiffness arising from the elastin-related parameters E and  and inversely proportional to  and , which define the maximum lung vital capacity through the distension of collagen. We know of no other quantitative model that addresses how changes in tissue constituents affect the shape of the PV curve.

However, we did not explicitly include the effects of surface tension in the model, and it may be that changes in surface tension in these diseases may contribute to changes in K. Surface tension contributes a significant portion to the elastic recoil of the lung clearly shown by the difference between PV curves for saline- and air-filled lungs (32, 35). Furthermore, Haber et al. (13) found that L_m correlated well with K in air-filled lungs but was unrelated to K determined from saline-filled lungs, suggesting that changes in air-space size were related to surface tension effects rather than tissue properties. However, we note as per the analysis of Hoppin and Hildebrandt (15) that surface tension properties are mechanically coupled. Removal of surface tension forces (e.g., by saline filling) alters alveolar wall geometry by changing tissue distention and thus affects the recoil attributable to the tissues.

The shapes of PV or length-tension curves can be altered by digestion of tissue constituents by proteases (16, 29). Karlinsky et al. (16) found that the chord compliance was increased in isolated hamster lungs exposed to elastase, but no change in vital capacity occurred. By contrast, collagenase caused an increase in vital capacity and a change in the shape of the PV curve. These observations are consistent with our model. However, when the measurements were done in vivo by the same authors, there was an increase in vital capacity observed in the elastase-treated lungs. This contradicts the theory that the collagen matrix is responsible for limiting the maximum extension. They attributed the difference between the in vitro and in vivo results to the repair process after elastase injury, which may disrupt the collagen matrix. Sata et al. (29) found that the stress-strain curves of microtissue specimens from hamster parenchyma could be described by two exponential functions, where the exponent of the first exponential was primarily responsible for the low-strain portion of the curve, whereas the upper portion of the curve was due to the exponent of the second exponential. They found that elastase treatment affected only the first exponent, whereas collagenase affected only the second component, which is in agreement with our model.

In any case, we do not claim that our model accounts for all aspects of the lung PV curve but merely that it embodies the important contribution from the tissues themselves whose significance for lung recoil is well accepted. Indeed, Wilson et al. (42) invoked our mechanism to explain the regional mechanical heterogeneity they observed during expirations from total lung capacity in dogs when using parenchymal markers. Such regional variation is predicted by our model, as different lung regions would have different samplings of the underlying collagen and elastin fiber distributions and thus would follow different local PV curves.

Concluding Comments

ACKNOWLEDGEMENTS

This work was supported by the Natural Sciences and Engineering Research Council of Canada, Fonds pour la Formation de Chercheurs et l'Aide à la Recherche (FCAR), the Medical Research Council of Canada, and the J. T. Costello Memorial Research Fund. J. H. T. Bates is a Chercheur-boursier of the Fonds de la Recherche en Santé du Québec. G. N. Maksym is supported by FCAR.

FOOTNOTES

Address for reprint requests: J. H. T. Bates, Meakins-Christie Laboratories, McGill Univ., 3626 St. Urbain St., Montreal, Quebec, Canada H2X 2P2 (E-mail: jason{at}meakins.lan.mcgill.ca).

Received 3 May 1996; accepted in final form 31 July 1996.

REFERENCES