|Home | About | Journals | Submit | Contact Us | Français|
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.
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.
The π-periodic von Mises distribution (Mardia and Jupp, 2000),
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)
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/π.
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,
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
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
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 θ.
The second integral is simpler and thus evaluated first. By making the substitution u = (θ-μ), the integral is converted to
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
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
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
Recognizing that a weighted sum of cosines and sines is equivalent to a single scaled and shifted cosine function, the integral is rewritten as
Incorporating the pre-integration terms and introducing
the final closed-form expression results:
And by similar methods, the other stress components can be written as
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.
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
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.
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.
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).
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.
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.
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.
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
The components derived from the expressions in the main body are
It is noted that only five components of D are independent.
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.