Abstract
We have developed a model of forces developed in lung tissue in which the stressbearing units are heterogeneous. Each element of the fiber network is composed of an idealized elastin and collagen element in parallel. Elastin is represented by linear springs and collagen by stiff strings that extend without resistance until taut. The model can quantitatively account for the nonlinear shape of the lengthtension curve of lung tissue strips when the knee lengths of the collagen fibers are distributed according to an inverse power law. The novel feature of this model is that as macroscopic strain increases the load is carried by progressively fewer elements with progressively higher forces, and preferential pathways of force transmission emerge within the matrix. The topology of these selforganizing pathways of force transmission takes the rough appearance of cracks, but, unlike real cracks, they represent the locus of force concentration rather than force release.
 lung tissue mechanics
 collagen
 elastin
 stressstrain
 power law
the elastic properties of lung tissue are determined by its stressbearing constituents, the collagen and elastin fibers. Since the seminal work by Mead et al. (21), it has been recognized that the static elastic properties of lung tissue are derived not only from the properties of the underlying tissue constituents but also from the way these constituents are organized. The collagen and elastin fibers are distributed in a meshlike network and vary in length, width, orientation, and curvature (19, 22, 25, 29). As the tissue is loaded, the amount of distension or strain arises from both stretching of individual fibers and their reorganization. It is commonly thought that elastin fibers elongate, whereas curved collagen fibers straighten and serve to stiffen the tissue and limit the final distension (10, 20, 28). Many models have been proposed for the mechanics of lung tissue incorporating the ideas of either fiber curvature, orientation or recruitment, or a combination of these effects. One common approach is to treat the tissue as a continuum where the mechanics are derived by integrating over the constituent properties (56, 1415, 17, 30, 32). Otherwise, a discrete model is constructed, with elements representing the combined roles of collagen and elastin as well as surface tension in simplified geometric configurations, usually representing a single alveolus or alveolar duct (3, 7, 9, 1113).
In this report, we have examined how forces develop in a twodimensional model for lung tissue elasticity with heterogeneously distributed elements. Each element of the model consists of an idealized collagen and elastin fiber in parallel, which extend and reorient in response to the applied uniaxial load. The collagen fibers are flaccid at low strain and only contribute to the elasticity of the elements when they become straight. In this study, we randomly assign the straightening lengths of the collagen fibers from various distributions and compare the forcelength curves computed to those measured in real lung tissue. We show that this leads to a heterogeneous distribution of forces among the elements and results in the development of preferential pathways of force transmission.
METHODS
The model.
We constructed the simulated tissue as a twodimensional finiteelement mesh consisting of interconnected line elements. Each line element represents a springstring pair. The springs represent elastin fibers, whereas the strings represent collagen fibers. The stressstrain curve for each line element is shown in Fig. 1. Each springstring pair begins at resting length with no load, and the strings are flaccid. As the tension in an individual pair increases, the unit extends according to the elastin stiffness (k
_{1}). At the knee length (l
_{k}), the string becomes taut, and the unit now extends according to the much stiffer parallel combination of the string and spring with combined stiffness (k
_{2}) (the string stiffness is thus k
_{2}−k
_{1}). Thus the lengthtension relationship of a fundamental unit is given by
Simulation. The tissuesimulation software was written in the Oberon programming language and environment (27). The solution was obtained by using the method of steepest descents, which makes use of the force gradient at the node positions. The solution algorithm used is as follows. The tension within each of the elements at each iteration is computed according toEq. 1
, allowing the resultant force on each node (intersection of elements) to be determined. The nodes are then moved down their respective force gradient according to
We simulated tissue strips with lengthtowidth ratios similar to those of the real lung tissue strips that we have reported previously (17). Models of 30 × 5 nodes (381 elements), 60 × 10 nodes (1,661 elements), and 90 × 15 nodes (3,841 elements) were simulated. The properties of the elements were chosen according to the following four schemes: 1) identicall
_{k} with fixedk
_{1} andk
_{2};2) identicall
_{k} with Gaussiandistributedk
_{1} andk
_{2}, keeping the same mean k
_{1} andk
_{2} as in1; 3) a powerlawdistributedl
_{k}, distributed between limits of 10^{−5} and 10^{5} according to
For Gaussiandistributed spring constants, no spring constant was permitted to be less than zero, andk _{2} was restricted to be greater thank _{1} within a single springstring pair. Because the units of length and force are arbitrary, length was normalized in all simulations to the resting length of the springs. These lengths were all identical. Spring constants were assigned stiffnesses in the range of 0.01 to 1,000 in units of force per unit of extension. The Gaussian random numbers used to assign values tok _{1 }andk _{2 }were generated according to the transformation method, and powerlawdistributed random numbers were chosen according to the rejection method (26). The node positions, element forces, and average axial force were monitored during the iterative process. Forcestrain curves were generated by stepping the tissue from a strain of zero in steps of 0.05 to a strain of 1.25 in 25 increments, allowing the tissue to converge at each step.
Calculation of stressstrain curves from the tissue simulations took from 3 min to 16 h on an IBM RS6000 computer, depending on the number of elements, distribution ofl _{k}, and degree of nonlinearity. In general, the greater the ratio ofk _{2} tok _{1}, the longer the convergence time. Tissues in which thel _{k} were distributed according to power laws posed severe numerical difficulties. This was primarily due to the fact that the element property distributions necessarily spanned several decades, and thus sampling the distributions adequately required a very large number of units. The largest simulations we attempted were 90 × 15 nodes. However, solution of some of these models required several days, particularly for high ratios ofk _{2} tok _{1}. Simulations of models having 60 × 10 nodes usually required less than a few hours, and so we used this size of model most often. Because the model was finite, different realizations chosen from the same element property distributions produced slightly different results. The amount of variation depended on b(Eq. 4 ), and was approximately inversely related to the number of elements. Nevertheless, the median forcelength curves from models of different sizes were very similar, so in the present study we report those simulations producing the median force at 0.60 strain.
RESULTS
When the model elements are all the same (l _{k} = 1.5,k _{1} = 0.1,k _{2} = 10), the tissue narrows rapidly at strains greater than ∼0.5, which corresponds to the level at which the elements begin to exceed their knee lengths, and elements near the ends of the tissue begin to carry more stress than central elements (Fig. 2). This behavior is due to the geometric distortion caused by the fixed endnode positions. Thus, although the element properties are all identical, the forces born by the elements are not (Figs. 2 and3) because of geometrical constraints.
When the spring constants (k _{1} andk _{2}) are chosen from Gaussian distributions (the mean values ofk _{1} andk _{2} are the same as before and have 50% SDs of their respective means as in Fig. 2), the forces carried by the elements vary considerably (Fig.4). However, despite the heterogeneous distribution of forces, the forcelength curve is nearly identical to that of tissue with homogeneously distributed elements (Fig.5). The forcestrain curves are essentially piecewise linear, with a “knee” at about strain 0.7. Thus neither simulation (Figs. 2 and 4) is representative of real lung tissue, which has a smoothly curvilinear lengthtension relationship. Instead, the simulated curves resemble those of a single element (Fig. 1), except that the knee is shifted rightward from strain 0.5 to 0.7. This is attributable to the additional strain afforded by changes in orientation, before the elements reach their knee lengths.
Figure 6 shows forcestrain curves from 60 × 10 node tissue strips withk _{1} = 0.1,k _{2} = 100, andl _{k} distributed according to Eq. 4 with values forb ranging from 0.9 to 1.1. This sensitivity was due to the rate at which the elements exceeded their respective l _{k} as the tissue was stretched and depended on both thel _{k} distribution and the value ofk _{1}. The value ofk _{2} was less important because, as an element reached its knee length, it was more likely to distort the nearby tissue rather than move up to the steep portion of its lengthtension curve. The forcelength curves were found to match the experimental data only when thel _{k} values were distributed with b(Eq. 4 ) very close to unity. Figure7 shows the forcelength curve from two realizations (b = 1), together with a stressstrain curve from a single dog lung tissue strip from Maksym and Bates (17).
The effect on the forcelength curve of distributing the spring constants according to Gaussian distributions as well as distributing the l _{k} according to a power law (b = 1) is shown in Fig. 8. Unlike the sensitivity to the exponent b of thel _{k} distribution, there was robustness to distributing bothk _{1} andk _{2} while keeping the l _{k}hyperbolically distributed. The lengthtension curve showed little change in shape even when the SDs of the stiffness distributions were 50% of their means (Fig. 8).
As the tissue (power lawl _{k},b = 1) is stretched, the degree of spatial heterogeneity is increased, and the degree of force heterogeneity is also increased (Fig. 9). At lower strains, elements with shorterl _{k} begin to carry progressively larger forces than do their neighbors, causing greater local strains. As the tissue lengthens, elements nearby these regions begin to form chains representing preferential pathways of force transmission. The corresponding forcelength curve calculated from this tissue strip is the lower of the two simulated curves in Fig. 7.
DISCUSSION
The key finding of this report is that in a network model with randomly assigned elemental properties there emerged selforganizing chains of force transmission. The model accounted for the lengthtension behavior of lung tissue using elements embodying only the essential mechanical characteristics of the collagen and elastin fibers, where the knee length of each element was randomly chosen according to an inverse powerlaw distribution.
Setnikar (28) first proposed the idea that elastin fibers were responsible for maintaining tension at lower strains while collagen took over the load at higher strains and limited the maximum distension (20). This is based on the reasoning that elastin is very extensible, rupturing at extensions exceeding 230% of its original length, whereas collagen can only extend 1–3% once taut and has a Young’s modulus 1,000 to 10,000 times larger than that of elastin. Thus deformation of lung tissue was implicitly ascribed to the stretching of elastin and the reorganization of collagen fibers. Each element in our model embodies this concept as a parallel pair of elements: an elastinlike spring and a collagenlike string, where the element “knee lengths” are distributed. The variation among elements invokes the natural variation of the structural properties of the fibers, which are known to be distributed over a wide range of lengths, widths, and curvatures (19, 22, 25, 29). Indeed, the inverse powerlaw distribution of knee lengths is similar in shape to the collagen fiber curvature distribution found for human alveolar wall, which shows a decreasing tail for the majority of the data (29). The deformation and reorganization of the fibers, as well as reorganization of large scale structures such as the alveolar wall and large septal boundaries under an applied load, lead to the bulk elasticity of lung tissue. However, the nature by which tensions are distributed throughout the tissue has not been well explored. We have developed a model that not only mimics the lengthtension behavior of lung tissue but also uses the underlying distributed nature of its constituents to examine force development within the tissue.
The nature by which forces develop in this model is our most significant finding (Fig. 9). As the tissue was stretched, some elements began to carry more load than did their neighbors. This led to local distortions that increased the strain on adjacent elements. This increased the likelihood that neighboring elements would become stressed to the stiff part of their lengthtension curve. With increasing stretch, small selforganizing chains of force transmission emerged and were distributed heterogeneously throughout the tissue. As the stretch continued, the chains began to link together, forming ever longer preferential pathways of force transmission. For large stretch, a continuous link formed from one end of the tissue to the other, where each element of the link was stressed beyond itsl _{k}. When thel _{k} were distributed hyperbolically, the progressive manner in which the chains formed and distributed tension led to forcelength curves mimicking actual lung tissue (Fig. 7). This may reflect the way in which collagen fibers are recruited with increasing strain to produce the smoothly stiffening stressstrain curve seen in real lung tissue.
The morphology and growth of these selforganizing force chains resemble that of crack propagation through a solid (18). However these force chains are really the opposite of cracks, or “anticracks.” Cracks form where regions of high stress cause the material to yield locally, thereby relieving the stress. By contrast, the anticracks in Fig. 9 represent regions where force is concentrated. Just as cracks propagate by increasing local distortions and stresses at the crack tip, these anticracks also grow by increasing local distortions, which then increase the load on the neighboring elements. The anticracks also selforganize by linking to nearby elements in a manner analogous to the selforganization shown in a recent model of fracture formation (4). Like the scale invariance exhibited in the roughness and topology of crack formation (33), lung tissue is similarly a rich source of scaleinvariant characteristics, within physiological limits. The lung shows powerlaw behavior in its branching network, in powerlaw avalanches of opening pressures (31), and in its impedance (2). It may be that the inverse powerlaw distribution (b = 1) of knee lengths found in our model is elicited due to a scaleinvariant character governing lung tissue organization and function.
Power laws are features of systems governed according to the theory of selforganized criticality (SOC), which has been applied to such phenomena as avalanches in sandpiles and earthquakes (1). In SOC, energy is continually fed into a system that is made of many connected energystoring elements, each with a maximum storage limit. When this limit in an element is reached, the energy spills over to the neighboring elements, bringing them closer to their respective energy limits. Periodically, a cascade (avalanche) of energy spills occurs. Eventually, the system reaches a critical state reflected in the spatial and temporal distributions of avalanches. These distributions are invariably described by power laws. Our tissue model also exhibits avalanchetype behavior in the sense that each elastic element stores energy up to a length limitl _{k}, at which point it suddenly becomes much stiffer. This means that additional strain must be absorbed mostly by neighboring elements while they are still shorter thanl _{k}. This leads to the formation of concentrations of force transmission, with force distributions (Fig. 3 C) that appear distributed as an inverse power law with exponent equal to 2. However, the similarities between our model and SOC systems end there, because our model is purely elastic with no mechanism for energy dissipation. Thus it cannot evolve into a critically stable state and is therefore not an SOC system. Nevertheless, our observations suggest that powerlaw distributions may be a general feature of systems that tend to selforganize, regardless of whether they are dissipative.
This formation of preferential pathways of force transmission is strikingly similar to how forces are transmitted through a sand bed. Liu et al. (16) studied stressinduced optical birefringence in a uniformly loaded cylinder containing identical pyrex spheres and found that stress is distributed exponentially and is organized into branching chains of force. These chains are oriented randomly, due to the randomness in the contact angles between beads, which is highly reminiscent of the patterns of the anticracks in our model that resulted from the randomness of the properties of model elements.
Is there evidence for a heterogeneous spatial or force distribution within the lung? The degree of spatial distortion in Fig. 9 is not normally observed, since most samples are fixed at defined pressures (25, 34). Although spatial distortion is observed in tissue sections exposed to contractile agonists (8, 24), this may arise from either heterogeneity in the lung constituents, differences in smooth muscle cell response, their distribution, or a combination of these effects. Inhomogeneity in the passive mechanics of lung tissue has been observed in vivo by using parenchymal markers by Wilson et al. (35), who showed that the regional pressurevolume curves follow markedly disparate paths. Conceivably, they were sampling different regions of the underlying distributions of the constituent properties.
Although spatial heterogeneity in constituent properties in the model led to increased distortion and a heterogeneous force distribution, this does not imply highly heterogeneous distributions of stress within the lung parenchyma or necessarily large amounts of tissue distortion. Higher forces can be maintained by thicker fibers, which would give a more even distribution of stresses. Indeed, the distributed nature of this model is founded on the fact that the lung tissue fiber properties are widely distributed (19, 25, 29). This raises the interesting possibility that a mechanism may exist for remodeling the connective tissue matrix in such a way that variations in stress are minimized. Mechanically transduced remodeling may be relevant to understanding lung tissue growth and maintenance and to diseases that compromise lung tissue elasticity, such as emphysema or fibrosis. Fibroblasts in vitro require external tension to be active and produce extracellular matrix proteins (23). The large local forces arising in our model are associated with large strains over short distances and are presumably a likely locus for fibroblast activation and matrix synthesis. Thus the preferential pathways of force transmission could lead to tissue growth exactly at the locus of the anticracks, thereby eventually relieving local stresses.
Acknowledgments
We thank Whitney deVries and William Thorpe for their assistance and their expert programming advice.
Footnotes

Address for reprint requests: J. H. T. Bates, MeakinsChristie Laboratories, McGill University, 3626 St. Urbain St., Montréal, Québec H2X 2P2, Canada (Email: jason{at}meakins.lan.mcgill.ca).

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, the Medical Research Council of Canada, the J. T. Costello Memorial Research Fund, and the National Institutes of Health.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.
 Copyright © 1998 the American Physiological Society