PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Biomech. Author manuscript; available in PMC 2010 July 22.
Published in final edited form as:
PMCID: PMC2707510
NIHMSID: NIHMS112903

A Closed-Form Structural Model of Planar Fibrous Tissue Mechanics

Abstract

Structural models of tissue mechanics, in which the tissue is represented as a sum or integral of fiber contributions for a distribution of fiber orientations, are a popular tool to represent the complex mechanical behavior of soft tissues. A significant practical challenge, however, is evaluation of the integral that defines the stress. Numerical integration is accurate but computationally demanding, posing an impediment to incorporation of structural models into large-scale finite element simulations. In this paper, a closed-form analytic evaluation of the integral is derived for fibers distributed according to a von Mises distribution and an exponential fiber stress-strain law.

Keywords: soft-tissue mechanics

Introduction

One of the major challenges in the mathematical description of soft-tissue mechanics is the lack of an obvious choice of constitutive equation that can capture the complex architecture and non-linear behavior of the tissue. An important conceptual advance, put forward by Lanir (Lanir, 1982; Lanir, 1983) and subsequently implemented and modified by many others (e.g., (Baek et al., 2007; Billiar and Sacks, 2000; Gasser et al., 2006; Pinsky et al., 2005; Sacks, 2003; Wagenseil and Okamoto, 2007)), is to introduce contributions from a continuous matrix and from a population of fibers, with the constituents assumed to deform affinely with the macroscopic deformation and the stresses assumed to be additive. Although those assumptions may not be valid in all circumstances (Chandran and Barocas, 2006; Hepworth et al., 2001; Sacks, 2003), the general form remains an attractive option and has been shown to capture a wide range of material behavior with relatively few fitting parameters. Structural information, in the form of fiber orientations, can be introduced based on optical methods (Sacks, 2003; Sander et al., to appear).

There are two basic strategies that have been employed to describe the fiber population. In the discrete approach, a few individual fiber populations, each aligned perfectly in a certain direction, are considered. In the continuous approach, a continuous distribution function (e.g., von Mises or a modified normal distribution) is used to represent the fiber population statistically. The continuous approach has the advantage that it applies to systems in which the fiber population is highly distributed with respect to direction, but the disadvantage that the total stress from the population must be expressed as an integral rather than as a sum of a small number of values. Evaluation of the integral, which is generally performed numerically, can be quite costly, encouraging an attempt to derive invariants that can capture the behavior without requiring the integral evaluation (Freed et al., 2005). Both that paper and (Gasser et al., 2006) derive analytical closures for a kinematic integral (i.e., an average measure of stretch in the fibers based on the fiber orientation distribution and the macroscopic kinematics), then apply a stress-strain formulation to the resulting quantity. In the case of (Freed et al., 2005), it is advocated that the strain energy should be expressed as the integral of the fiber stress-strain function from unit stretch (i.e., unloaded) to the averaged stretch. In (Gasser et al., 2006), an additional invariant is added to the strain energy function to incorporate the extra kinematic parameter. Neither approach involves direction integration of the stress over the distribution of fiber directions.

In this paper, a closed-form analytical expression based on the exact integral of the stress in a distribution of fibers is derived for an exponential fiber constitutive equations and a von Mises distribution. Other constitutive equations and fiber distribution equations could be considered, but the choices are made because of both have been used in the past and because they combine to produce a particularly simple result.

Model Construction

Fiber Distribution

The π-periodic von Mises distribution (Mardia and Jupp, 2000),

f(θ;κ,μ)=1πI0(κ)exp{κcos[2(θμ)]},θ[0,π)
(1)

where κ is a concentration parameter, μ is the preferred direction, and I0 is the modified Bessel function of the first kind and order 0, is a good choice to describe a planar collection of fibers because (1) it is naturally restricted to the range [0,π), (2) it is a limiting distribution in certain cases (Mardia and Jupp, 2000), and (3) it allows direct calculation of certain moments by means of the identity (Abramowitz and Stegun, 1964)

In(z)=1π0πezcosθcos(nθ)dθ=1π0πezcos2θcos(2nθ)dθ
(2)

which provides the normalization factor in (1). As κ becomes large, f approaches a normal distribution with variance 1/(4κ), and as κ becomes small, f approaches a uniform distribution on the half circle, f = 1/π.

Fiber Constitutive Equation

A critical feature of the structural approach is the transfer of constitutive description from the tissue level to the fiber level. That is, the model requires a constitutive equation for the fiber rather than for the tissue. An exponential form based on the Green strain of the fiber is a popular choice (Billiar and Sacks, 2000; Driessen et al., 2007) and used here as follows,

Sf=A{exp[B(λf21)]1}
(3)

where A and B are constitutive constants, Sf is the second Piola-Kirchoff stress in the fiber, and λf is the fiber stretch for a fiber aligned with angle θ, defined kinematically by

λf2=C1cos2(θα)+C2sin2(θα)=C11+cos[2(θα)]2+C21cos[2(θα)]2=C1+C22+C1C22cos[2(θα)]
(4)

where C1 and C2 are the squared principal stretches, i.e., the eigenvalues of the Green deformation tensor C = FTF, and α is the angle corresponding to the direction of the first eigenvector. Although other possibilities could be considered, the constitutive model of equation (3) was chosen because of its popularity and because of the simple closed-form integral that can be derived. Although the constants A and B do not have direct physical meaning, they capture the overall stiffness (A) and nonlinearity (B) of the fiber response; that is, increasing A makes the fiber stiffer at all strains, and increasing B makes the fiber response more sharply nonlinear. It is further noted that in the limit of small strain (λf → 1), equation (3) reduces to a linear response with modulus 2AB.

Since the stress is the sum of all fiber stresses, and the fibers are assumed to be stretched affinely with the macroscopic deformation, the total second Piola-Kirchoff stress Sij is given by

Sij=0πei(θ)ej(θ)Sf(θ;C1,C2,α)f(θ;κ,μ)dθ
(5)

where ei(θ) is the i component of the unit vector in the fiber direction θ, Sf is the fiber stress defined by equation (3), and f is the von Mises distribution. Although the coordinate directions remain arbitrary at this point in the analysis, it is noted that e1(θ) = cos θ and that e2(θ) = sin θ.

Analytical Evaluation of the Integral

Substituting Sf from equation (3) into (5), we write

S11=AπI0(κ)0πcos2θ{exp[B(λf21)]1}exp{κcos[2(θμ)]}dθ=AπI0(κ)0πcos2θexp{B(λf21)+κcos[2(θμ)]}dθAπI0(κ)0πcos2θexp{κcos[2(θμ)]}dθ
(6)

The second integral is simpler and thus evaluated first. By making the substitution u = (θ-μ), the integral is converted to

μπμcos2(u+μ)exp{κcos[2u]}dθ=0πcos2(u+μ)exp{κcos[2u]}du
(7)

where the change in the limits of integration is possible because the intregrand is π-periodic. Expressing cos2(u+μ) in terms of cos[2(u+μ)] and then applying the cosine-of-a-sum formula gives

0π1+cos2μcos2usin2μsin2u2exp{κcos[2u]}du=120πexp{κcos[2u]}du+cos2μ20πcos2uexp{κcos[2u]}dusin2μ20πsin2uexp{κcos[2u]}du=π2I0(κ)+πcos2μ2I1(κ)+0=π2[I0(κ)+I1(κ)cos2μ]
(8)

where the first two integrals are evaluated using (2), and the third is zero because the integrand is symmetric about u = π/2.

The first integral from (6) is evaluated as follows. Expanding λf2 from equation (4) and applying the cosine-of-a-difference formula, the integral to be evaluated (neglecting the pre-factor temporarily) is expressed as

0πcos2θexp{B(C1+C22+C1C22cos[2(θα)]1)+κ(cos2θcos2μ+sin2θsin2μ)}dθ=exp(BC1+C221)0πcos2θexp{B(C1C22cos[2(θα)])+κ(cos2θcos2μ+sin2θsin2μ)}dθ
(9)

where the second expression is obtained by moving terms independent of θ outside the integral. When cos[2(θ-α)] is expanded and the pre-integration factors are lumped together, the integral can be rewritten concisely as

0πcos2θexp{k1cos2θ+k2sin2θ}dθ
(10)

where

k1B2(C1C2)cos2α+κcos2μk2B2(C1C2)sin2α+κsin2μ
(11)

Recognizing that a weighted sum of cosines and sines is equivalent to a single scaled and shifted cosine function, the integral is rewritten as

0πcos2θexp{γcos2(θδ)}dθ
(12)

where

γk12+k22δ12tan1(k2k1)
(13)

By the same argument as applied in (7)-(8) above, the integral in (12) is

π2[I0(γ)+I1(γ)cos2δ]
(14)

Incorporating the pre-integration terms and introducing

k0A2exp[B(C1+C221)]
(15)

the final closed-form expression results:

S11=k0I0(γ)I0(κ)+k0I1(γ)I0(κ)cos2δA2AI1(κ)2I0(κ)cos2μ
(16)

And by similar methods, the other stress components can be written as

S22=k0I0(γ)I0(κ)k0I1(γ)I0(κ)cos2δA2+AI1(κ)2I0(κ)cos2μ
(17)
S12=k0I1(γ)I0(κ)sin2δAI1(κ)2I0(κ)sin2μ
(18)

Note that when C1 = C2 = 1 (no extension), k0 = A, γ = κ, and δ = μ, so the stress is, as required, zero. Differentiation of the above expressions is obviously cumbersome but not difficult, making use of the chain rule and the derivative relations for Bessel's functions (and observing that A, B, κ, and μ are independent of deformation). The tangent stiffness matrix is provided in the Appendix.

Determination of Structural Parameters

If an entire probability density function is available (e.g., by light scattering (Bowes et al., 1999; Sacks, 2003) or by Fourier analysis of an image (Marquez, 2006; Sander and Barocas, 2008)), then the von Mises distribution can be regressed to the data to determine κ and μ, or multiple distributions can be added if the data are multimodal. Other methods, such as polarized-light microscopy (Tower et al., 2002), provide only a measure of degree of anisotropy, which, in the simplest approximation (Chandran and Barocas, 2003), is related to the quantity

0πcos[2(θμ)]f(θ)dθ
(19)

or I1(κ)/I0(κ) for the von Mises distribution. Although the inverse Bessel functions are not generally available, the function I1(κ)/I0(κ) is monotonic in κ, so numerical estimation of κ from the degree of anisotropy data is straightforward, especially since the approximation I1(κ)/I0(κ) ≈ κ/2 provides a very good initial guess for κ < 1.

Application to Planar Biaxial Extension of a Bioartificial Tissue

Although the underlying structural model has been applied many times and thus should not need further evaluation, a simple test is conducted to confirm the validity of the closed-form expressions and to demonstrate their application. Jhun et al. (accepted) recently performed displacement-controlled equibiaxial strain experiments on the homogeneous central region of compacted fibroblast-populated collagen cruciforms. The cruciforms had different arm widths, with greater difference in arm width leading to greater anisotropy. In this paper, results for arm width ratios of 1:0.5 and 1:0.375 were analyzed by fitting three parameters in equations (16) and (17) for each arm width ratio: κ, A, and B. The fourth parameter, μ, was set to 0 since the loading direction corresponded to the characteristic direction of collagen fiber alignment, as determined by quantitative polarized light microscopy. The Cauchy stress tij was calculated from the second Piola-Kirchoff stress by assuming incompressibility of the sample and no shear (α = 0), so t11 = λ12S11 and t22 = λ22S22.

Figure 1 shows the experimental data of Jhun et al. and least-squares fits of the model. As expected given the past success of such models at predicting mechanical behavior in equibiaxial loading, the fit is excellent. The final model parameters were κ = 0.742, A = 372 kPa, and B = 3.57 for the 1:0.375 case (Figure 1a), and κ = 0.444, A = 542 kPa, and B = 3.01 for the 1:0.5 case (Figure 1b). Although it is reiterated that the purpose of this paper is to demonstrate the analytical solution, not to test the model itself, certain observations are in order. First, one would expect that A and B would be nearly the same for the two cases if the underlying collagen populations are the same except for their orientation; differences can be attributed to sample variation, different degrees of compaction between the two cases, sample inhomogeneity, the role of cells in the mechanics of the sample, or any other confounding factor that is not captured in the model. Second, Jhun reported the birefringence of the sample to be 33.3±0.97°/mm for the 1:0.375 case and 18.9±0.72°/mm for the 1:0.5 case; the ratio of the two is 1.76±0.08. Calculating I1(κ)/I0(κ) for the two fits and taking the ratio yields a value of 1.61, slightly lower but not unreasonable given the multiple assumptions made in deriving the model.

Figure 1
Fit of Model to Data of Jhun et al. (2008)

Simulation of Uniaxial Extension of an Obliquely Aligned Sample

As an example of implementation of the model in a finite-element setting, we have modeled uniaxial extension of a sample with the fibers aligned obliquely to the pull direction (A = 540 kPa, B = 3, κ = 0.6, μ = 60°). Figure 2 shows the deformed sample with the principle stretches (Figure 2a) and principal Cauchy stresses (Figure 2b). The results show that the off-axis alignment of the fibers lead to skewing of the strain field. The problem was solved twice, once using equations (16)-(18) for the stresses, and once using numerical integration of equation (5). The Gnu Scientific Library was used for both the Bessel functions and the quadrature (Gauss-Kronrod quadrature; on average, 43 points were sufficient for an error of 10-4). The results were the same in both cases, but the Bessel function was about 200 times faster for the entire calculation (i.e., evaluation of the integrals as well as other processes associated with the calculation).

Figure 2
Simulated Uniaxial Stretch of an Obliquely Aligned Sample

Discussion

The closed-form solution derived above provides the potential for more rapid solution of computational mechanics problems associated with fiber-population-based structural models. We found a speedup of about 200, but caution should be taken because results will depend on the specific methods used, especially for the Bessel functions. There is a wide range of approaches to the complex problem of describing tissue mechanics, and the utility of this derivation is restricted to cases where the underlying model is appropriate. Even within the class of structural models, certain fiber constitutive forms may not allow analytical integration. It is also important to recognize that the von Mises distribution does not describe all data sets, so other options might need to be considered. Nevertheless, the commonness of the von Mises distribution and the potential for summation of multiple distributions is a strong point in favor of this approach.

The angle δ and magnitude γ represent the coupling between the fiber affine stretch and fiber population distributions. Specifically, the angle α, which depends solely on the macroscopic deformation and not on the fiber distribution, is the direction in which a fiber would be stretched the most by a given macroscopic deformation. The angle μ, in contrast, depends solely on the fiber distribution and not on the macroscopic deformation; it is the direction in which a fiber is most likely to be oriented. The angle δ depends on both the deformation and the distribution, and it is the direction that contributes most to the total stress (more precisely, to the portion of the stress arising from the first integral in equation (6)). That such a quantity emerged naturally from the derivation is physically satisfying. The magnitude γ, as can be seen by comparing equations (1) and (12), functions similarly to the von Mises concentration parameter κ. That is, a large γ indicates that most of the stress is coming from fibers in a narrow orientation band, either because that is where most of the fibers are or because that is where they are most stretched (or, likely, a combination of the two); a small γ indicates that fibers in many directions are all contributing to the total stress. As a demonstration, the orientation, fiber stress, and γ-δ distributions are shown in Figure 3 for the structural parameters fitted to the 1:0.375 data and an imposed 10% : 1% biaxial stretch with the large stretch at an angle 45° from the principal fiber direction. The angle δ is 12.3° - at larger angles, even though the fibers are stretched more, there are fewer of them, and at smaller angles, the many fibers are not stretched as much.

Figure 3
Angular Distributions in Structural Model

In comparing the current model with the closed-form models of others (Freed et al., 2005; Gasser et al., 2006), the major difference is that the current model provides a macroscopic constitutive equation (16)-(18) based entirely on the fiber mechanics and without any additional constitutive assumptions at the tissue scale. While the direct link from fiber to tissue scale is intuitively appealing, it is not clear that the model is actually better than the others – all approaches require a considerable amount of simplification of the complex mechanical behavior of a tissue, and all reduce to a reasonable set of equations and fitting parameters. The choice of specific model strategy is best made by considering the ability of any given model to capture the mechanics of the tissue of interest.

One important drawback of the current model in contrast to others is its restriction to two dimensions. While it may be possible to work out analytical solutions for certain spherical distributions, the current model is applicable only to planar tissues. Curved membranes could be modeled by the appropriate coordinate transformation, but a fully three-dimensional distribution of fibers is not modeled by this method. Likewise, there is no accounting for incompressibility of the tissue, either within the plane or through the plane. One could account for incompressibility by adding in an incompressible isotropic matrix to the stress (Driessen et al., 2007; Stylianopoulos and Barocas, 2007) or to the strain energy (Gasser et al., 2006).

Finally, it is observed that equations (16)-(18) could be viewed not as a fiber model per se but as a specialized constitutive model. The equations could be fit to experimental data using κ, μ, A, and B as fitting parameters as was done in the example (possibly estimating κ and μ from other experiments) without making a specific claim that the equation represents an underlying fiber population. Such an approach would be similar to the Arruda-Boyce model and its variations (Bischoff et al., 2000; Kuhl et al., 2005; Palmer and Boyce, 2008), which use a simple unit-cell representation of fiber mechanics – the model, although not a perfect description of the micromechanical environment, is able to describe a wide range of material responses with a relatively small number of parameters. The parameters in the current model represent four different effects: κ is the degree of anisotropy, μ is the preferred direction, A is the stiffness of the fiber, and B is the nonlinearity of the fiber stress-strain response.

Supplementary Material

01

Acknowledgments

The author thanks Choon-Sik Jhun for providing the data analyzed in this paper and Prof. Henryk Stolarkski for assistance in the implementation of the code. This work was supported by the National Institutes of Health (R01EB005813, R01HL071538) and by a resources grant from the Minnesota Supercomputing Institute.

Appendix

For finite element implementation, use of Newton-Raphson iteration to solve the nonlinear stress balance requires an expression for the tangent stiffness (i.e., the derivative of stress with respect to strain). Although we omit the details of the derivation for the sake of brevity, we provide here the tangent stiffness matrix (detailed derivation is available in supplementary material). Taking Sij to be the 2nd Piola-Kirchhoff stress defined in the main article, and Ekl to be the Green-Lagrange strain, Ekl = (Ckl - δkl), the tangent stiffness is defined by

DijklSijEkl

The components derived from the expressions in the main body are

D1111=k0B2I0(κ)[I0(γ)(1+cos22δ)+I1(γ)(2cos2δcos4δγ)]
D1122=D2211=12D1212=k0B2I0(κ)[I0(γ)(sin22δ)+I1(γ)(cos4δγ)]
D1112=2D1211=k0B2I0(κ)[I0(γ)(sin4δ)+I1(γ)(2sin2δ2sin4δγ)]
D2212=2D1222=k0B2I0(κ)[I0(γ)(sin4δ)+I1(γ)(2sin2δ+2sin4δγ)]
D2222=k0B2I0(κ)[I0(γ)(1+cos22δ)+I1(γ)(2cos2δcos4δγ)]

It is noted that only five components of D are independent.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  • Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. U.S. Govt. Printing Office; 1964. p. 1046.
  • Baek S, Valentin A, Humphrey JD. Biochemomechanics of cerebral vasospasm and its resolution: II. Constitutive relations and model simulations. Annals of Biomedical Engineering. 2007;9:1498–1509. [PubMed]
  • Billiar KL, Sacks MS. Biaxial mechanical properties of the native and glutaraldehyde- treated aortic valve cusp: Part II - A structural constitutive model. Journal of Biomechanical Engineering-Transactions of the Asme. 2000;4:327–335. [PubMed]
  • Bischoff JE, Arruda EM, Grosh K. Finite element modeling of human skin using an isotropic, nonlinear elastic constitutive model. Journal of Biomechanics. 2000;6:645–652. [PubMed]
  • Bowes LE, Jimenez MC, Hiester ED, Sacks MS, Brahmatewari J, Mertz P, Eaglstein WH. Collagen fiber orientation as quantified by small angle light scattering in wounds treated with transforming growth factor-beta(2) and its neutalizing antibody. Wound Repair and Regeneration. 1999;3:179–186. [PubMed]
  • Chandran PL, Barocas VH. Affine vs. non-affine fibril kinematics in collagen networks: theoretical studies of network behavior. Journal of Biomechanical Engineering. 2006:259–270. [PubMed]
  • Chandran PL, Barocas VH. Microstructural mechanics of collagen gels in confined compression: poroelasticity, viscoelasticity, and collapse. Journal of Biomechanical Engineering. 2003;2:152–166. [PubMed]
  • Driessen NJ, Mol A, Bouten CV, Baaijens FP. Modeling the mechanics of tissue-engineered human heart valve leaflets. Journal of Biomechanics. 2007;2:325–334. [PubMed]
  • Freed AD, Einstein DR, Vesely I. Invariant formulation for dispersed transverse isotropy in aortic heart valves: an efficient means for modeling fiber splay. Biomechanics and Modeling in Mechanobiology. 2005;(23):100–117. [PubMed]
  • Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of the Royal Society, Interface / the Royal Society. 2006;6:15–35. [PMC free article] [PubMed]
  • Hepworth DG, Steven-fountain A, Bruce DM, Vincent JFV. Affine versus non-affine deformation in soft biological tissues, measured by the reorientation and stretching of collagen fibres through the thickness of compressed porcine skin. Journal of Biomechanics. 2001;3:341–346. [PubMed]
  • Jhun CS, Evans MC, Barocas VH, Tranquillo RT. Fiber reorientation in planar biaxial loading of bioartificial tissues possessing prescribed alignment. Journal of Biomechanical Engineering accepted. [PMC free article] [PubMed]
  • Kuhl E, Garikipati K, Arruda EM, Grosh K. Remodeling of biological tissue: Mechanically induced reorientation of a transversely isotropic chain network. Journal of the Mechanics and Physics of Solids. 2005;7:1552–1573.
  • Lanir Y. Constitutive equations for fibrous connective tissues. Journal of Biomechanics. 1983;1:1–12. [PubMed]
  • Lanir Y. Constitutive equations for fibrous connective tissues. Journal of Biomechanics. 1982;1:1. [PubMed]
  • Mardia KV, Jupp PE. Directional Statistics. Wiley; 2000. p. 429.
  • Marquez JP. Fourier analysis and automated measurement of cell and fiber angular orientation distributions. International Journal of Solids and Structures. 2006:6413–6423.
  • Palmer JS, Boyce MC. Constitutive modeling of the stress-strain behavior of F-actin filament networks. Acta Biomaterialia. 2008;3:597–612. [PubMed]
  • Pinsky PM, van der Heide D, Chernyak D. Computational modeling of mechanical anisotropy in the cornea and sclera. Journal of Cataract and Refractive Surgery. 2005;1:136–145. [PubMed]
  • Sacks MS. Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar-collagenous tissues. Journal of Biomechanical Engineering-Transactions of the Asme. 2003;2:280–287. [PubMed]
  • Sander EA, Tranquillo RT, Barocas VH. Image-Based Biomechanics of Collagen-Based Tissue Equivalents: Multiscale Models Compared To Fiber Alignment Predicted by Polarimetric Imaging. IEEE Eng in Medicine and Biology to appear.
  • Sander EA, Barocas VH. Comparison of 2D fiber network orientation measurement methods. Journal of Biomedical Materials Research. Part A 2008 [PubMed]
  • Stylianopoulos T, Barocas VH. Multiscale, structure-based modeling for the elastic mechanical behavior of arterial walls. Journal of Biomechanical Engineering. 2007;4:611–618. [PubMed]
  • Tower TT, Neidert MR, Tranquillo RT. Fiber alignment imaging during mechanical testing of soft tissues. Annals of Biomedical Engineering. 2002;10:1221–1233. [PubMed]
  • Wagenseil JE, Okamoto RJ. Modeling cell and matrix anisotropy in fibroblast populated collagen vessels. Biomechanics and Modeling in Mechanobiology. 2007;3:151–162. [PubMed]