Soft composite based hyperelastic model for anisotropic tissue characterization

Abstract

Soft tissues are complex anisotropic composite systems comprising of multiple differently oriented layers of fiber embedded within a soft matrix. To date, soft tissues have been mainly characterized using simplified linear elastic material models, isotropic viscoelastic and hyperelastic models, and transversely isotropic models. In such models, the effect of fiber volume fraction (FVF), fiber orientation, and fiber-matrix interactions are missing, inhibiting accurate characterization of anisotropic tissue properties. The current work addresses this literature gap with the development of a novel soft composite based material framework to model tissue anisotropy. In this model, the fiber and matrix are considered as separate hyperelastic materials, and fiber-matrix interaction is modeled using multiplicative decomposition of the deformation gradient. The effect of the individual contribution of the fibers and matrix are introduced into the numerical framework for a single soft composite layer, and fiber orientation effects are incorporated into the strain energy functions. Also, strain energy formulations are developed for multiple soft composite layers with varying fiber orientations and contributions, describing the biomechanical behavior of an entire anisotropic tissue block. Stress-strain relationships were derived from the strain energy equations for a uniaxial mechanical test condition. To validate the model parameters, experimental models of soft composites tested under uniaxial tension were characterized using the novel anisotropic hyperelastic model (R² = 0.983). To date, such a robust anisotropic hyperelastic composite framework has not been developed, which would be indispensable for experimental characterization of tissues and for improving the fidelity of computational biological models in future.

Keywords

Tissue anisotropy characterization composite hyperelastic

Introduction

Soft tissues such as the skin, skeletal muscles, connective tissues, and tissues forming organs (such as the brain and myocardial tissues) are typically inhomogeneous^1,2 and exhibit regional and directional anisotropy in three dimensions.^3–5 Figure 1 shows a schematic of the various components of a soft tissue. Tissue anisotropy is mainly attributed to the variations in distribution of collagen fibers² from one site of the body to another. In the past, collagen fiber distribution in human cadavers and animal models were studied using dissection and histological investigations.^6,7 Advancement in the imaging techniques in the recent years have allowed to look at fiber distributions in the human body using Diffusion Tensor Magnetic Resonance Imaging (DT-MRI) and confocal imaging.⁸ In literature, there have been very few attempts to model tissue anisotropy in Finite Element (FE) models involving organs.^1,3,9–14 Most other FE and experimental models so far have used isotropic and transversely isotropic hyperelastic models for tissue material characterization.

Figure 1.

Schematic of the various components of a soft tissue.

The material behavior of an isotropic soft material (Figure 2(a)) is characterized using the strain-energy function ( $ψ$ ), which depends on the type of material.¹⁵ In an isotropic hyperelastic model, the strain-energy function is dependent on the principal stretches ( $λ_{1}$ , $λ_{2}$ , and $λ_{3}$ ) or the Cauchy-green tensor invariants ( $I_{1}$ , $I_{2}$ and $I_{3}$ , which are also functions of the principal stretches)¹⁵ as shown in equations (1) and (2). Some of the well-known isotropic hyperelastic constitutive models used to characterize soft tissues are Neo-Hookean, Mooney-Rivlin, Yeoh, Ogden, Humphrey, Martins and Veronda-Westmann.¹⁵

ψ_{Isotropic} = ψ (I_{1}, I_{2}, I_{3})

(1)

I_{1} = \sum_{i = 1}^{3} λ_{i}^{2}, I_{2} = \sum_{i, j = 1}^{3} λ_{i}^{2} λ_{j}^{2}, I_{3} = \prod_{i = 1}^{3} λ_{i}^{2}

(2)

Figure 2.

Schematic of soft composite stretching in: (a) Isotropic hyperelastic model and (b) transversely-Isotropic hyperelastic model.

Mooney-Rivlin model is the first hyperelastic material model developed, which was successful in precise polynomial curve fitting of mechanical stress-strain test data obtained for non-linear isotropic rubber-like materials. Equation (3) presents this model, where the applied stress in a uniaxial test is in terms of ‘ $λ$ ’, the uniaxial stretch and c₁ and c₂, the two material constants which are predicted using curve fitting. Mooney-Rivlin model has been used extensively in the past to predict soft tissue material behavior under uniaxial and biaxial testing.^16–25 Yeoh model, shown in equation (4) is a reduced polynomial form of Mooney-Rivlin model, which can accurately predict the non-linear behavior of incompressible rubber-like materials. Yeoh model has been used recently to predict the mechanical properties and brain and breast tissues.^26,27

σ_{Mooney} = 2 (λ^{2} - \frac{1}{λ}) (c_{1} + c_{2} \frac{1}{λ})

(3)

σ_{Yeoh} = 2 (λ^{2} - \frac{1}{λ}) (c_{1} + 2 c_{2} (I_{1} - 3) + 3 c_{3} {(I_{1} - 3)}^{2})

(4)

Neo-Hookean model (equation (5)), is a primitive form of hyperelastic model similar to the Hooke’s law, developed in the year 1948 to study the mechanical properties of vulcanized rubber. It has been observed to be accurate for strains less than 20% and fail in predicting material behavior under biaxial stresses. Neo-Hookean model has been superseded by Mooney-Rivlin model in accuracy of prediction, however it is still been used to predict non-linear soft tissue behavior.²⁸ Ogden model (equation (6)) predicts the mechanical properties for incompressible, isotropic and stress rate independent rubber-like materials more accurately than Mooney-Rivlin model and is also stable at high strains, and has been used often to characterize biological tissues.^20,29–34

σ_{Hookean} = 2 (λ^{2} - \frac{1}{λ}) c_{1}

(5)

\begin{array}{l} σ_{Ogden} = c_{1} (λ^{c_{2}} - 2^{- 1 + c_{2}} λ^{- c_{2} / 2}) \\ + c_{3} (λ^{c_{4}} - 2^{- 1 + c_{4}} λ^{- c_{4} / 2}) \\ + c_{5} (λ^{c_{6}} - 2^{- 1 + c_{6}} λ^{- c_{6} / 2}) \end{array}

(6)

Humphrey and Yin in 1987^35,36 developed a composite based material model for myocardial soft tissues. The soft tissue was assumed to be made of various families of non-interacting fibers and a homogeneous matrix, both of which were assumed to undergo gross deformation. A pseudo-strain energy function was proposed which depends on limited structural information, and both on the matrix and fiber contribution information, which were obtained through multiaxial stress-strain experiments. Humphrey and Yin’s model (equation (7)), which is an exponential function with two material constants, was able to accurately predict anisotropic, incompressible and non-linear behavior of soft tissues.^37–42

σ_{Humphrey} = 2 (λ^{2} - \frac{1}{λ}) c_{1} c_{2} e^{c_{2} (I_{1} - 3)}

(7)

Martin’s hyperelastic model^15,43 was introduced recently as a modified form of the transversely isotropic (Figure 2(b)) and incompressible hyperelastic model by Humphrey and Yin. An extra term signifying the passive strain energy stored in the muscle fibers due to the uniaxial experiments, was added to the overall strain energy expression. The resulting stress-strain relationship given by equation (8) was able to predict soft tissue material behavior, and was subsequently used by Martins et al. to predict the non-linear behavior of pelvic tissues.^43–50 Veronda-Westmann model⁵¹ (equation (9)) is another composite based model like Humphrey and Yin’s model, in which the energy density function depends on the contribution of the fibers and the matrix material. This model was developed based on experimental testing on cat’s skin, and has been used to predict the material behavior of liver, breast tissues and silicon rubber phantoms.^27,52–55

\begin{array}{l} σ_{Martins} = 2 (λ^{2} - \frac{1}{λ}) c_{1} c_{2} e^{c_{1} (I_{1} - 3)} + 2 λ (λ - 1) c_{3} c_{4} e^{c_{3} {(λ - 1)}^{2}} \end{array}

(8)

σ_{Veronda - Westmann} = 2 (λ^{2} - \frac{1}{λ}) c_{1} c_{2} (e^{c_{2} (I_{1} - 3)} - \frac{1}{2 λ})

(9)

To date, the isotropic hyperelastic models are able to model the gross mechanical properties of a soft tissue, and does not consider the effect of fibers. The transversely-isotropic hyperelastic material models consider the effect of fibers, and in some cases the fiber contribution. However, most of such models lack the consideration of fiber-matrix interactions, fiber orientations and distributions in layers. In this work, a novel soft composite based hyperelastic material model was developed which considers a soft tissue as an arrangement of several differently oriented fiber family layers within a matrix material. This model considers the individual fiber and matrix contributions in each layer, the effect of varying fiber orientations across layers, and the fiber-matrix interactions. The next section discusses the soft composite hyperelastic model development. Then, the validation of the model against the mechanical properties of soft composite materials developed in literature is presented, which is followed by conclusions in the final section.

Anisotropic hyperelastic model

Numerical model

To include the anisotropic effect of fiber inclusions in a matrix (soft composite model) with all fibers organized initially in the direction of vector $a_{0}$ (Figure 3), the energy function takes the form $ψ = ψ ({C,a}_{0})$ , where $C$ is the Lagrangian right Cauchy-Green deformation tensor. To use the invariant approach for modeling transverse isotropy, two additional invariants related to $a_{0}$ , are defined as $I_{4} = a_{0} .C {.a}_{0} = λ_{F}^{2}, I_{5} = a_{0} {.C}^{2} {.a}_{0}$ where $λ_{F}$ is the fiber stretch.^42,56–59 The strain energy function is rewritten as $ψ (C, a_{0}) = ψ (I_{1}, I_{2}, I_{3}, I_{4}, I_{5})$ . To simplify computations, both the matrix and the fibers are assumed to be incompressible materials. As the stiffness of the matrix and fibers are different, so deformation will be non-uniform in matrix and fibers. We use deformation tensor $C_{m}$ in matrix to define the invariants in strain energy function of the matrix given by $ψ^{M} = C_{10}^{m} (I_{1}^{m} - 3)$ where $I_{1}^{m} = t r (C_{m})$ and the material model is neo-Hookean type.^42,56–59 $C_{10}^{m}$ is related to shear modulus of the matrix material. The strain energy stored in the fibers is modeled as a function of $I_{4}$ , which is based on the assumption that the strain energy stored in the fibers depends mainly on the fiber elongation, where $λ_{F} = \sqrt{I_{4}}$ is the actual fiber stretch. To account for non-linear dependence of strain energy on fiber stretch, and also to account for the three-dimensional aspect of the fiber, the strain energy stored in fibers is written as $ψ^{F} = C_{10}^{m} f (I_{4}) (I 1 f - 3)$ where $I 1 f = t r (C_{f})$ and the material model is neo-Hookean type, $C_{f}$ is the deformation tensor of fibers and $C_{10}^{m} f (I_{4})$ is related to stretch-dependent shear modulus of the fiber.^42,56–59

Figure 3.

(a) Schematic of fiber rotation due to soft composite deformation and (b) overall soft composite stretching due to fiber-matrix interaction, in the novel soft composite based anisotropic hyperelastic model formulation

Considering initial fiber direction $a_{0}$ along $x_{1}$ direction of rectangular Cartesian coordinate system (Figure 3), i.e. $a_{0} = e_{1} = {[\begin{array}{l} 1 & 0 & 0 \end{array}]}^{T}$ . Assuming second arbitrary perpendicular axis $e_{2}$ such that $e_{1} {.e}_{2} = 0$ . Third axis $e_{3} = e_{1} \otimes e_{2} = 0$ . With a deformation $F$ , the deformed fiber direction ${Fa}_{0}$ is rotated back by angle $γ$ to original fiber direction $a_{0} = e_{1}$ . Let the rotation tensors for this be $R_{1}$ which is independent of choice of $e_{2}$ , and $R_{2}$ which depends on the $e_{2}$ axis. Using $e_{2}$ as the rotation axis, $R_{1} {Fe}_{2}$ is rotated back to $e_{1} {-e}_{2}$ plane. The new deformation gradient tensor $F^{*}$ is written as $F^{*} {= R}_{2} R_{1} F$ . As $R_{1}$ and $R_{2}$ are rigid body rotations, ${(F^{*})}^{T} F^{*}$ is the same as $F^{T} F$ . $F^{*}$ is used as starting point of the following multiplicative decomposition of deformation gradient.^42,56–59 From definition of $F^{*}$ in the current deformation considered in Figure 3, $F_{21}^{*} = F_{31}^{*} = F_{32}^{*} = 0$ and $F_{11}^{*} = | {Fa}_{0} | = λ_{F}$ . $F^{*}$ in $e_{1} {,e}_{2} {,e}_{3}$ coordinate system is given by equation (10).^42,56–59

F^{*} = [\begin{array}{l} λ_{F} & F_{12}^{*} & F_{13}^{*} \\ 0 & F_{22}^{*} & F_{23}^{*} \\ 0 & 0 & F_{33}^{*} \end{array}]

(10)

The deformation gradient $F^{*}$ is multiplicatively decomposed into uniaxial deformation gradient $F_{f}^{*}$ in the fiber direction and subsequent shear deformation gradient $F_{s}^{*}$ such that $F^{*} {= F}_{s}^{*} F_{f}^{*}$ . $F_{f}^{*}$ and $F_{s}^{*}$ are given by equation (11).^42,56–59

\begin{array}{l} F_{f}^{*} = diag [\begin{array}{l} λ_{F} & λ_{F}^{- \frac{1}{2}} & λ_{F}^{- \frac{1}{2}} \end{array}], \\ F_{s}^{*} = F^{*} {(F_{f}^{*})}^{- 1} = [\begin{array}{l} 1 & F_{12}^{*} λ_{F}^{\frac{1}{2}} & F_{13}^{*} λ_{F}^{\frac{1}{2}} \\ 0 & F_{22}^{*} λ_{F}^{\frac{1}{2}} & F_{23}^{*} λ_{F}^{\frac{1}{2}} \\ 0 & 0 & F_{33}^{*} λ_{F}^{\frac{1}{2}} \end{array}] \end{array}

(11)

By suitable choice of $e_{2}$ axis, $F_{22}^{*} λ_{F}^{\frac{1}{2}} = 1$ . $F_{12}^{*}$ and $F_{13}^{*}$ are along fiber shears, and $F_{23}^{*}$ is termed as transverse shear. Shear deformations are related to invariants of $C$ given by equation (12).^42,56–59

{(F_{12}^{*})}^{2} + {(F_{13}^{*})}^{2} = \frac{I_{5} - I_{4}^{2}}{I_{4}}, {(F_{23}^{*})}^{2} = I_{1} - \frac{I_{5} + 2 \sqrt{I_{4}}}{I_{4}}

(12)

Modeling the effect of fiber and matrix contributions

Soft tissues have been studied extensively in literature and have been found to exhibit a wide range of stress-strain behavior due to differences in collagen fiber content and fiber distribution.^15,52,57 The highest collagen fiber content have been observed in connective tissues making them significantly stronger than other tissues.⁶⁰ Additionally, due to diseases such as Pelvic Organ Prolapse (POP), soft (pelvic) tissue changes have been attributed to collagen fiber variations.³ To model the effect of individual fiber and matrix contribution on a soft composite material property, the matrix volume fraction (MVF) $v_{m}$ and fiber volume fraction (FVF) $v_{f}$ terms were introduced. The strain energy stored in the whole composite during uniaxial deformation $F_{f}^{*}$ is given by equation (13),^42,56–59 where $I_{1} (F_{f}^{*}) = λ_{F}^{2} + \frac{2}{λ_{F}} = I_{4} + 2 I_{4}^{- \frac{1}{2}}$ is the first invariant of $C_{f}^{*} = {(F_{f}^{*})}^{T} F_{f}^{*}$ .

\begin{array}{l} ψ_{f} = ψ_{f}^{M} + ψ_{f}^{F} = v_{m} C_{10}^{m} (I_{1} (F_{f}^{*}) - 3) \\ + v_{f} C_{10}^{m} f (I_{4}) (I_{1} (F_{f}^{*}) - 3) \end{array}

(13)

On application of shear deformation $F_{s}^{*}$ , fiber stretch and stiffness remains constant. In this step, the composite is treated as neo-Hookean matrix embedded with neo-Hookean fibers where $f (I_{4})$ is ratio of the fiber and matrix stiffnesses. The quantity of shear deformation is computed by adding the terms as shown in equation (14), where $I_{1} (F_{f}^{*}) = I_{4} + 2 I_{4}^{- \frac{1}{2}}$ .^42,56–59

\begin{array}{l} {(F_{12}^{*})}^{2} + {(F_{13}^{*})}^{2} + {(F_{23}^{*})}^{2} = \frac{I_{5} - I_{4}}{I_{4}} + I_{1} - \frac{I_{5} + 2 \sqrt{I_{4}}}{I_{4}} \\ = I_{1} - I_{1} (F_{f}^{*}) \end{array}

(14)

The energy stored during deformation $F_{s}^{*}$ of the neo-Hookean composite is given by equation (15), where $c_{3}$ is the effective neo-Hookean modulus, and $I_{1} - I_{1} (F_{f}^{*})$ is the portion of $I_{1}$ remaining after deformation .^42,56–59

ψ_{s} = c_{3} [I_{1} - I_{1} (F_{f}^{*})]

(15)

Now strain energy function for the anisotropic composite model is sum of the strain energies of the two steps (equations (13) and (15)), which is given by equation (16).^42,56–59 In equation (16), $C_{10}^{m}$ and $C_{10}^{m} f (I_{4})$ are substituted with constants $c_{1}$ and $c_{2}$ respectively.

\begin{array}{l} ψ = ψ_{f} + ψ_{s} = v_{f} c_{2} [I_{1} (F_{f}^{*}) - 3] + v_{m} c_{1} [I_{1} (F_{f}^{*}) - 3] \\ + c_{3} [I_{1} - I_{1} (F_{f}^{*})] \end{array}

(16)

For Humphrey model fibers and matrix definitions, the modified strain energy functions and the net strain energy functions are given by equations (17) and (18) respectively.¹⁵

ψ^{M} = c_{1} [e^{c_{2} (I_{1}^{m} - 3)} - 1], ψ^{F} = c_{3} [e^{c_{4} (I_{1}^{f} - 3)} - 1]

(17)

\begin{array}{l} ψ = v_{m} c_{1} [e^{c_{2} (I_{1} (F_{f}^{*}) - 3)} - 1] + v_{f} c_{3} [e^{c_{4} (I_{1} (F_{f}^{*}) - 3)} - 1] \\ + c_{5} [I_{1} - I_{1} (F_{f}^{*})] \end{array}

(18)

In equations (16) and (18), the term $[I_{1} - I_{1} (F_{f}^{*})]$ goes to zero with $I_{1} (F_{f}^{*}) = I_{1}$ as $λ_{F}$ (the stretch in the direction of the fiber) in case of a uniaxial tensile test coincides with the applied stretch $λ$ ^{15,42,56–59} for a fiber direction along the direction of stretch. Also, substituting $I_{1} (F_{f}^{*}) = λ^{2} + \frac{2}{λ}$ , equations (16) and (18) can be rewritten as equations (19) and (20) respectively.

ψ_{Neo - Hookean / Modified} = (v_{m} c_{1} + v_{f} c_{2}) (λ^{2} + \frac{2}{λ} - 3)

(19)

\begin{array}{l} ψ_{Humphrey / Modified} = v_{m} c_{1} [e^{c_{2} (λ^{2} + \frac{2}{λ} - 3)} - 1] \\ + v_{f} c_{3} [e^{c_{4} (λ^{2} + \frac{2}{λ} - 3)} - 1] \end{array}

(20)

The principal Cauchy stress for a given hyperelastic strain energy function was found using equation (21), where $σ_{2} = σ_{3} = 0$ are the boundary conditions for the uniaxial tension experiment.^{15,42,56–59} Here $λ_{1} = λ$ and $λ_{3} = 0$ .

σ_{1} = λ_{1} \frac{\partial ψ}{\partial λ_{1}} - λ_{3} \frac{\partial ψ}{\partial λ_{3}}, σ_{2} = σ_{3} = 0

(21)

Substituting the strain energy functions from equations (19) and (20) into equation (21), the stress versus stretch equations are given by equations (22) and (23).

σ_{Neo - Hookean / Modified} = 2 (v_{m} c_{1} + v_{f} c_{2}) (λ - \frac{1}{λ^{2}})

(22)

\begin{array}{l} σ_{Humphrey / Modified} = 2 v_{m} c_{1} c_{2} (λ - \frac{1}{λ^{2}}) [e^{c_{2} (λ^{2} + \frac{2}{λ} - 3)}] \\ + 2 v_{f} c_{3} c_{4} (λ - \frac{1}{λ^{2}}) [e^{c_{4} (λ^{2} + \frac{2}{λ} - 3)}] \end{array}

(23)

Material characterization is conducted using equations (22) and (23) by substituting the fiber contribution ( $v_{f}$ ) and matrix contribution ( $v_{m}$ ) in the stress-stretch equations and estimating the values of the constants c’s.

Modeling the effect of fiber orientation

To account for single layer fiber orientation effect in the soft composites (Figure 4), $λ_{F} = λ \cos θ$ was substituted in equation (16), so that $I_{1} (F_{f}^{*}) = λ_{F}^{2} + \frac{2}{λ_{F}} = λ^{2} \cos^{2} θ + \frac{2}{λ \cos θ}$ , $[I_{1} - I_{1} (F_{f}^{*})] = [λ^{2} (\sin^{2} θ) + \frac{2}{λ} (1 - \sec θ)]$ and $0 \leq θ < 90^{o}$ . It should be noted that the $θ$ representation used here for the fiber angle is different from the $γ$ used in the ‘Numerical model’ section to explain fiber rotation. The modified stress-strain relationships of types neo-Hookean (equation (24)) and Humphrey (equation (25)) were derived similar to equations (17) to (23), starting from equation (16), which includes the effect of fiber angle $θ$ , fiber deformation, and fiber rotation $γ$ .

\begin{array}{l} σ_{Neo - Hookean / Modified 2} = 2 (v_{m} c_{1} + v_{f} c_{2}) \\ \times (λ \cos θ - \frac{1}{λ^{2} \cos^{2} θ}) \\ + 2 c_{3} [λ (\sin^{2} θ) - \frac{1}{λ^{2}} (1 - \sec θ)] \end{array}

(24)

\begin{array}{l} σ_{Humphrey / Modified 2} = 2 v_{m} c_{1} c_{2} (λ \cos θ - \frac{1}{λ^{2} \cos^{2} θ}) \\ \times [e^{c_{2} (λ^{2} \cos^{2} θ + \frac{2}{λ \cos θ} - 3)}] \\ + 2 v_{f} c_{3} c_{4} (λ \cos θ - \frac{1}{λ^{2} \cos^{2} θ}) \\ \times [e^{c_{4} (λ^{2} \cos^{2} θ + \frac{2}{λ \cos θ} - 3)}] \\ + 2 c_{5} [λ (\sin^{2} θ) - \frac{1}{λ^{2}} (1 - \sec θ)] \end{array}

(25)

Figure 4.

Schematic showing the stretch along the fiber ( $λ_{F}$ ) which is at an angle $+ θ$ or $- θ$ with respect to the soft composite uniaxial stretch $λ$ .

Figure 5.

Schematic for novel anisotropic material formulation for soft composites with a single matrix and multiple fiber families.

Modeling the effect of multiple fiber layers at arbitrary orientations

For a soft composite with a single matrix and multiple fiber families embedded within it (Figure 5), the mechanical property depends only on the cumulative effect of fiber directions.^42,56–59 For a composite with $n$ fiber families where the fibers are oriented at angles $θ_{1}, θ_{2}, θ_{3}, \dots θ_{i} \dots θ_{n}$ (where $0 \leq θ_{i} < 90 °$ and $n > 1$ ), volume fractions are $v_{m}$ for the matrix and $v_{f 1}, v_{f 2}, v_{f 3}, \dots v_{f i} \dots v_{f n}$ for each fiber family respectively, the homogenized strain energy function can be written as equation (26).^42,56–59

ψ_{Soft Composite} = \sum_{i = 1}^{n} ψ^{i}

(26)

The modified generic stress-strain relationships based on equations (24) to (26) for a soft tissue composite with multiple fiber families are given by equations (27) and (28).

\begin{array}{l} σ_{Soft Composite : Neo - Hookean} \\ = \sum {\begin{cases} 2 ((1 - v_{f i}) c_{1} + v_{f i} c_{2}) (λ \cos θ_{i} - \frac{1}{λ^{2} \cos^{2} θ_{i}}) \\ + 2 c_{3} [λ (\sin^{2} θ_{i}) - \frac{1}{λ^{2}} (1 - \sec θ_{i})] \end{cases}}, \\ 1 \leq i \end{array}

(27)

\begin{array}{l} σ_{Soft Composite : Humphrey} \\ = \sum {\begin{cases} 2 (1 - v_{f i}) c_{1} c_{2} (λ \cos θ_{i} - \frac{1}{λ^{2} \cos^{2} θ_{i}}) [e^{c_{2} (λ \cos θ_{i} - \frac{1}{λ^{2} \cos^{2} θ_{i}} - 3)}] \\ + 2 v_{f i} c_{3} c_{4} (λ \cos θ_{i} - \frac{1}{λ^{2} \cos^{2} θ_{i}}) [e^{c_{4} (λ^{2} \cos^{2} θ_{i} + \frac{2}{λ \cos θ_{i}} - 3)}] \\ + 2 c_{5} [λ (\sin^{2} θ_{i}) - \frac{1}{λ^{2}} (1 - \sec θ_{i})] \end{cases}}, \\ 1 \leq i \leq n \end{array}

(28)

Results and model validation

For validation of the current model developed for characterization of the material properties of tissues, soft composite materials in literature⁶¹ were employed. These soft composite materials have been observed to closely mimic the anisotropic material behavior of brain tissue,^62,63 pelvic tissues,^64,65 and the skin.^66–68 The true stress versus true stretch observed with soft composites with FVF of 0.17, 0.35, 0.52, 0.61, and 0.78 respectively, were traced and plotted in Figure 6. For each plot, the fiber contribution ( $v_{f}$ ) and matrix contribution ( $v_{m}$ ), and single fiber angle ( $θ$ ) of 0°, and number of layers ( $n = 1$ ) were substituted into equations (27) and (28), and $σ - λ$ curve fitting was conducted using the optimization tool in Microsoft Excel. The model parameters estimated using equations (27) and (28) are presented in Table 1. An average correlation index (R² value) of all the curve fits was estimated to be 0.98.

Figure 6.

Average stress versus stretch plots traced from the literature on mechanical testing of soft composites with varying FVF.⁶¹

Table 1.

Hyperelastic curve fit coefficients for the soft composite materials⁶¹ with varying FVF’s (based on Figure 6).

Hyperelastic curve fit coefficients
Fiber volume fraction (FVF)	Anisotropic soft composite (Neo-Hookean)		Anisotropic soft composite (Humphrey)
Fiber volume fraction (FVF)	c₁	c₂	c₁	c₂	c₃	c₄
0.17	0.95	1.23	0.53	0.41	0.40	0.61
0.35	1.80	2.18	1.15	0.33	1.20	0.63
0.52	2.65	2.80	0.72	0.44	1.22	0.64
0.61	3.24	3.56	0.73	0.42	1.97	0.65
0.78	3.96	4.15	0.81	0.45	2.31	0.65

A second set of soft composite properties were considered from literature,⁶¹ comprising of the ones with multiple fibers and varying fiber angles. Figure 7 shows the trace of the average true stress versus true stretch profiles of the uniaxially tested soft composites with fiber angles $θ = 0 °, \pm 45 °, 90 °$ and varying $v_{m}$ and $v_{f}$ (generated due to different number of fibers).

Figure 7.

Average stress versus stretch plots traced from the literature on mechanical testing of soft composites with one, two and three fiber configurations and different fiber angles.⁶¹

Table 2 presents the model parameters calculated for the second set of soft composites using the novel formulations (equations (27) and (28)). Table 2 also includes the FVF’s corresponding to each of the specimen configurations mentioned in Figure 7. The soft composite variants with specifically ±90° were avoided, as mentioned previously in the material modeling section, due to multiple terms in the constitutive model going to infinity. The soft composite with a 90° fiber angle in the current uniaxial test framework is assumed to behave exactly like a composite with FVF = 0. However, it should be mentioned that in reality, such a consideration is not true, as loadings are not uniaxial but multiaxial in nature, which will be subject of future studies. The average R² correlation index estimated for this set was 0.986. Using the average correlation index values from the hyperelastic modeling for the two set of soft composite experiments, the correlation index of the newly developed hyperelastic model was concluded to be 0.983.

Table 2.

Hyperelastic curve fit coefficients for the soft composite materials⁶¹ with different number of fibers and orientations presented in Figure 7.

Hyperelastic curve fit coefficients
Fiber (Orientation)	FVF	Anisotropic soft composite (Neo-Hookean)			Anisotropic soft composite (Humphrey)
Fiber (Orientation)	FVF	c₁	c₂	c₃	c₁	c₂	c₃	c₄	c₅
One Fiber (0)	0.21	0.05	3.23	–	0.31	0.55	0.75	0.75	–
Two Fiber (0)	0.38	0.01	2.95	–	0.45	0.65	0.85	0.65	–
Three Fiber (0)	0.56	0.01	2.37	–	0.55	0.75	1.14	0.61	–
One Fiber (±45)	0.16	–1.16	9.25	0.95	–1.15	4.85	4.87	5.13	1.3
Two Fiber (±45)	0.27	–2.54	9.41	1.33	–1.05	6.14	2.81	6.18	1.3
Three Fiber (±45)	0.41	–5.88	10.63	1.8	–2.13	5.59	2.52	5.84	2.6

Discussion

A review of the state of the art of isotropic and transversely-isotropic hyperelastic models in literature indicated the lack of capabilities to date for accurate modeling of soft tissue anisotropy. In the current work, a novel soft composite based hyperelastic characterization framework was developed to model tissue anisotropy due to different layers of fibers and fiber orientations within a tissue. This material model considered both the matrix and fibers as hyperelastic materials, and took into account the matrix contribution and the individual fiber contributions, the fiber-matrix interaction, and fiber orientations for each fiber family within the soft composite. Using this novel anisotropic soft tissue characterization model, stress-strain relationships were generated for uniaxial mechanical testing. The model formulations were validated with experiments on biofidelic soft composites, which have widely been used in literature to simulate accurate soft tissue properties.

There are several limitations of the novel anisotropic hyperelastic model. First, a fiber family was assumed to not transition or interact with the adjacent fiber families. In a natural soft tissue composite system, fibers run across the layers, the effect of which will be considered in future modeling studies. Second, within a layer, the fiber orientations also change from one location to another in a natural tissue. In the current model, this location-specific change in the orientation of the fiber families was also not considered. Histology and microscopic imaging needs to be conducted to study such fiber orientation changes in different tissue systems. Third, the mechanical effect of the presence of blood vessels, capillaries and fat within tissue layers was neglected. Inclusion of the fat layer specifically in future modeling studies may lead to better tissue simulation in the study of the effect of impact loads. Fourth, effects of strain rates on the response of soft composite layers were also not modeled. Fifth, similar to other hyperelastic material models, this model did not account for energy dissipation, and is not able to simulate viscoelastic properties such as hysteresis, relaxation, and creep, observed in soft tissues. Energy dissipation will be considered in future to further improve the current model. In order to realistically characterize the soft tissue material behavior using the current model under cyclic loading and unloading, the mechanical response can be modeled as hyperelastic with different material constants at loading and unloading respectively. This method suggested by Fung et al.⁶⁹ is known as pseudoelastic modeling, which utilizes the elastic theory to model an inelastic material, such as a soft tissue. Sixth and lastly, the novel computational model was tested mainly for uniaxial mechanical test conditions. In reality, soft tissue composite loads are biaxial and multi-axial, the effect of which will be considered as a part of future work.

To date, the soft composite hyperelastic model developed in this work most comprehensively captures tissue anisotropy. This model can improve experimental modeling of tissues, and better characterization of normal and diseased tissues at different sites on the body. The computational biomechanical modeling community will also greatly advantage from using this novel hyperelastic model for better tissue material modeling, and consequently improved accuracy of results from analyses. This robust model also opens up opportunities for multiscale modeling of tissues including the study of the effects of microstructural changes on macro scale tissue properties.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Arnab Chanda

References

Picinbono

Delingette

Ayache

Nonlinear and anisotropic elastic soft tissue models for medical simulation. In: IEEE international conference on robotics and automation, 2001. Proceedings 2001 ICRA. Piscataway: IEEE, 2001.

Holzapfel

GA.

Biomechanics of soft tissue. Handbook Mater Behav Models 2001; 3: 1049–1063.

Chanda

Unnikrishnan

Roy

, et al. Computational modeling of the female pelvic support structures and organs to understand the mechanism of pelvic organ prolapse: a review. Appl Mech Rev 2015; 67: 040801.

Chanda

Ghoneim

Pumping potential of a two-layer left-ventricle-like flexible-matrix-composite structure. Compos Struct 2015; 122: 570–575.

Ghoneim

Chanda

Pumping potential of a left-ventricle-like flexible matrix-composite structure. In: e-Proceedings of the 19th international conference on composite materials (ICCM19), Montreal, Canada, 28 July 2013.

Lowry

Gilligan

Katersky

EM.

The determination of collagen and elastin in tissues, with results obtained in various normal tissues from different species. J Biolog Chem 1941; 139: 795–804.

Neuman

Logan

MA.

The determination of collagen and elastin in tissues. J Biol Chem 1950; 186: 549–556.

Basser

Pajevic

Pierpaoli

, et al. In vivo fiber tractography using DT‐MRI data. Magn Reson Med 2000; 44: 625–632.

Kruger

Nash

, et al. Anisotropic effects of the levator ani muscle during childbirth. Biomech Model Mechanobiol 2011; 10: 485–494.

10.

Wang

Haynor

Kim

An investigation of the importance of myocardial anisotropy in finite-element modeling of the heart: methodology and application to the estimation of defibrillation efficacy. IEEE Trans Biomed Eng 2001; 48: 1377–1389.

11.

Franzone

Guerri

Pennacchio

, et al. Spread of excitation in 3-D models of the anisotropic cardiac tissue. II. Effects of fiber architecture and ventricular geometry. Math Biosci 1998; 147: 131–171.

12.

Carlsen

Daphalapurkar

NP.

The importance of structural anisotropy in computational models of traumatic brain injury. Front Neurol 2015; 6

13.

Annaidh

, et al. Characterization of the anisotropic mechanical properties of excised human skin. J Mech Behav Biomed Mater 2012; 5: 139–148.

14.

Chanda

Graeter

Unnikrishnan

Effect of blasts on subject-specific computational models of skin and bone sections at various locations on the human body. Springfield, MO: AIMS Press, 2015.

15.

Martins

Jorge

Ferreira

A comparative study of several material models for prediction of hyperelastic properties: application to silicone‐rubber and soft tissues. Strain 2006; 42: 135–147.

16.

Beatty

MF.

Topics in finite elasticity: hyperelasticity of rubber, elastomers, and biological tissues – with examples. Appl Mech Rev 1987; 40: 1699–1734.

17.

Chabanas

, et al. Comparison of linear and non-linear soft tissue models with post-operative CT scan in maxillofacial surgery. Medical simulation 2004: 19–27.

18.

Guerin

Elliott

DM.

Quantifying the contributions of structure to annulus fibrosus mechanical function using a nonlinear, anisotropic, hyperelastic model. J Orthop Res 2007; 25: 508–516.

19.

Hirokawa

Tsuruno

Three-dimensional deformation and stress distribution in an analytical/computational model of the anterior cruciate ligament. J Biomech 2000; 33: 1069–1077.

20.

Desai

JP.

Characterization of soft-tissue material properties: large deformation analysis. Medical simulation 2004: 28–37.

21.

Kim

Srinivasan

MA.

Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse FE parameter estimation. In: Medical Image Computing and Computer-assisted Intervention – MICCAI, Palm Springs, CA, USA, 26-29 October 2005, pp.599-606. Berlin: Springer. Berlin: Springer, 2005, pp.599–606.

22.

Lally

Dolan

Prendergast

Cardiovascular stent design and vessel stresses: a finite element analysis. J Biomech 2005; 38: 1574–1581.

23.

Lally

Reid

Prendergast

Elastic behavior of porcine coronary artery tissue under uniaxial and equibiaxial tension. Ann Biomed Eng 2004; 32: 1355–1364.

24.

Mendis

Stalnaker

Advani

A constitutive relationship for large deformation finite element modeling of brain tissue. J Biomech Eng 1995; 117: 279–285.

25.

Prendergast

Lally

Daly

, et al. Analysis of prolapse in cardiovascular stents: a constitutive equation for vascular tissue and finite-element modelling. J Biomech Eng 2003; 125: 692–699.

26.

Kaster

Sack

Samani

Measurement of the hyperelastic properties of ex vivo brain tissue slices. J Biomech 2011; 44: 1158–1163.

27.

O'Hagan

Samani

Measurement of the hyperelastic properties of 44 pathological ex vivo breast tissue samples. Phys Med Biol 2009; 54: 2557–2569.

28.

Miller

Method of testing very soft biological tissues in compression. J Biomech 2005; 38: 153–158.

29.

Holzapfel

Ogden

RW.

Biomechanics of soft tissue in cardiovascular systems. Vol. 441. Berlin: Springer, 2003.

30.

Velardi

Fraternali

Angelillo

Anisotropic constitutive equations and experimental tensile behavior of brain tissue. Biomech Model Mechanobiol 2006; 5: 53–61.

31.

El Sayed

Mota

Fraternali

, et al. A variational constitutive model for soft biological tissues. J Biomech 2008; 41: 1458–1466.

32.

Gao

Lister

Desai

JP.

Constitutive modeling of liver tissue: experiment and theory. Ann Biomed Eng 2010; 38: 505–516.

33.

Holzapfel

Gasser

Ogden

RW.

A new constitutive framework for arterial wall mechanics and a comparative study of material models. J Elast Phys Sci Solids 2000; 61: 1–48.

34.

Holzapfel

Gasser

Ogden

RW.

Comparison of a multi-layer structural model for arterial walls with a Fung-type model, and issues of material stability. J Biomech Eng 2004; 126: 264–275.

35.

Humphrey

Yin

On constitutive relations and finite deformations of passive cardiac tissue: I. A pseudostrain-energy function. J Biomech Eng 1987; 109: 298–304.

36.

Humphrey

Yin

A new constitutive formulation for characterizing the mechanical behavior of soft tissues. Biophys J 1987; 52: 563–570.

37.

May-Newman

Yin

A constitutive law for mitral valve tissue. J Biomech Eng 1998; 120: 38–47.

38.

Raghavan

Vorp

DA.

Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability. J Biomech 2000; 33: 475–482.

39.

Seshaiyer

Humphrey

JD.

A sub-domain inverse finite element characterization of hyperelastic membranes including soft tissues. J Biomech Eng 2003; 125: 363–371.

40.

Van Loocke

Lyons

Simms

A validated model of passive muscle in compression. J Biomech 2006; 39: 2999–3009.

41.

Wang

DHJ

Makaroun

Webster

, et al. Mechanical properties and microstructure of intraluminal thrombus from abdominal aortic aneurysm. J Biomech Eng 2001; 123: 536–539.

42.

Weiss

Maker

Govindjee

Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput Meth Appl Mech Eng 1996; 135 : 107–128.

43.

Martins

, et al. Modelling the mechanical behavior of soft tissues using hyperelastic constitutive models. In: II International Conference on Computational Bioengineering – ICCB, Universidade de Brasília, Brasília, Brazil, 15 July 2005.

44.

Calvo

Peña

Martins

, et al. On modelling damage process in vaginal tissue. J Biomech 2009; 42: 642–651.

45.

Martins

Peña

Calvo

, et al. Prediction of nonlinear elastic behaviour of vaginal tissue: experimental results and model formulation. Comput Meth Biomech Biomed Eng 2010; 13: 327–337.

46.

Martins

Peña

Jorge

RMN

, et al. Mechanical characterization and constitutive modelling of the damage process in rectus sheath. J Mech Behav Biomed Mater 2012; 8: 111–122.

47.

Martins

PALS

Jorge

RMN

Ferreia

AJM

, et al. Vaginal tissue properties versus increased intra-abdominal pressure: a preliminary biomechanical study. Gynecol Obstet Invest 2011; 71: 145–150.

48.

Martins

PALS

Filho

ALS

Fonseca

AMRM

, et al. Uniaxial mechanical behavior of the human female bladder. Int Urogynecol J 2011; 22: 991–995.

49.

Peña

Calvo

Martínez

, et al. Experimental study and constitutive modeling of the viscoelastic mechanical properties of the human prolapsed vaginal tissue. Biomech Model Mechanobiol 2010; 9: 35–44.

50.

Peña

Martins

Mascarenhas

, et al. Mechanical characterization of the softening behavior of human vaginal tissue. J Mech Behav Biomed Mater 2011; 4: 275–283.

51.

Veronda

Westmann

Mechanical characterization of skin – finite deformations. J Biomech 1970; 3: 111–124.

52.

Groves

Coulman

Birchall

, et al. An anisotropic, hyperelastic model for skin: experimental measurements, finite element modelling and identification of parameters for human and murine skin. J Mech Behav Biomed Mater 2013; 18: 167–180.

53.

Mehrabian

Campbell

Samani

A constrained reconstruction technique of hyperelasticity parameters for breast cancer assessment. Phys Med Biol 2010; 55: 7489–7508.

54.

Pavan

Madsen

Frank

, et al. Nonlinear elastic behavior of phantom materials for elastography. Phys Med Biol 2010; 55: 2679–2692.

55.

Yin

Sun

Wang

, et al. Modeling of elastic modulus evolution of cirrhotic human liver. IEEE Trans Biomed Eng 2004; 51: 1854–1857.

56.

Guo

Peng

Moran

A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus. J Mech Phys Solids 2006; 54: 1952–1971.

57.

Holzapfel

Ogden

RW.

Constitutive modelling of arteries. In: Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 31 March 2010, paper no. 2118, pp.1551–1597.

58.

Brown

Smith

LM.

A simple transversely isotropic hyperelastic constitutive model suitable for finite element analysis of fiber reinforced elastomers. J Eng Mater Technol 2011; 133: 021021.

59.

Holzapfel

GA.

Nonlinear solid mechanics. Vol. 24. Chichester: Wiley, 2000.

60.

Fullmer

HM.

A comparative histochemical study of elastic, pre-elastic and oxytalan connective tissue fibers. J Histochem Cytochem 1960; 8: 290–295.

61.

Chanda

Callaway

Tissue anisotropy modeling using soft composite materials. Appl Bionics Biomech 2018; 2018: 1–9.

62.

Chanda

, et al. Biofidelic human brain tissue surrogates. MAMS J 2016: 1335–1341.

63.

Chanda

Unnikrishnan

Human tissue simulants for study of traumatic brain injury (TBI). In: Proceedings of the American Society for Composites: Thirty-First Technical Conference, Williamsburg, VA, 19 September 2016.

64.

Chanda

Flynn

Unnikrishnan

Biomechanical characterization of normal and prolapsed vaginal tissue surrogates. J Mech Med Biol 2018; 18: 1750100.

65.

Chanda

Unnikrishnan

Flynn

Biofidelic vaginal tissue surrogate. US Patent App. 62/263,942, 2015.

66.

Chanda

Unnikrishnan

Flynn

Biofidelic human skin simulant. US Patent App. 62/189, 2015.

67.

Chanda

Unnikrishnan

Flynn

Biofidelic skin simulant. 2017. Google Patents.

68.

Chanda

Unnikrishnan

Flynn

, et al. Experimental study on tissue phantoms to understand the effect of injury and suturing on human skin mechanical properties. Proc Inst Mech Eng H 2017; 231: 80–91.

69.

Fung

YC.

Biomechanics: mechanical properties of living tissues. Berlin: Springer Science & Business Media, 2013.