Soft biological tissues (SBTs) are recognized as anisotropic hyperelastic materials since they are naturally made of fiber reinforcements and a fluid-like matrix for supporting reversible finite deformations. Rubber-like materials are often treated as isotropic hyperelastic materials. The normalized nominal stress-stretch curves in uniaxial tension for SBTs in a fiber direction and rubber-like materials are shown in Figure 1. In uniaxial tension tests of SBTs, less stiffening at initial
Figure 1. Normalized stress-stretch curves of soft tissue and rubber in uniaxial tension tests.
stretches and more stiffening at final stretches have been observed, and the stress-stretch curve behaves more linearly at both initial and final stretches, making SBTs far more efficient as compared to rubber-like materials. Anisotropic constitutive modeling for finite deformations of SBTs requires experimental characterizations combined with theoretical predictions. As a common practice, theoretical models are fitted with experimental data tested in certain deformation modes and the fitted models predict deformations in untested modes. Thus, the theoretical development of constitutive models, the optimal design of experimental tests, and their numerical implementations into finite element methods are crucial in analyses and designs associated with SBTs.
Kinematics, conservation principles, and constitutive relations are the three pillars of continuum mechanics. Fundamental continuum theories for modeling SBTs as anisotropic hyperelastic materials were presented in the monographs by Fung (1993)  , Holzapfel (2000)  , and Cowin (2013)  . In constitutive relations, the second Piola-Kirchhoff stress tensor is generally related to the right Cauchy-Green tensor through a stored energy density functional or simply a stored energy functional, which needs to be rigorously determined. Thus, the determination of an anisotropic CSE functional is indispensable for constitutive modeling of anisotropic hyperelastic materials, including SBTs and fiber reinforced polymers.
A CSE functional is the key to both anisotropic and isotropic constitutive modelings for hyperelastic materials. For anisotropic constitutive modelings, strain invariant approaches dominate the established anisotropic models for SBTs since material properties must be invariant from all angles of observation. This objectivity ensures frame-indifference. Polynomial, power, exponential, and other functions of invariant components have been used to construct ICGs. Anisotropic constitutive models for SBTs have been reviewed through different perspectives by Humphrey (2003)  , Gasser, Ogden, and Holzapfel (2006)  , Holzapfel and Ogden (2010)  , Chagnon, Rebouah, and Favier (2015)  , and Lanir (2017)  .
Finite deformations can be measured through the right Cauchy-Green tensor 
where is the deformation gradient tensor. The three invariants of right Cauchy-Green tensor, , , and , are defined by
where is the second-order unit tensor and the matrix operators, “:”, “tr”, “adj”, and “det”, denote double contraction, trace, adjugate, and determinant operations, respectively.
The adjugate, cofactor, and determinant of are related by the following equations
where the superscript T denotes the transpose operator of a matrix.
The first-order derivatives of invariants are given by
The second-order derivatives of invariants are given by
where is the fourth-order zero tensor. The derivatives of second-order symmetric tensors are defined as
in which , for example, is the Kronecker delta and
Anisotropic finite deformations can be modeled through invariant components of the right Cauchy-Green tensor and structural tensors. According to Spencer (1971, 1984)   , one invariant component in a preferred direction i, , can be defined by
where is a structural tensor, is a unit vector along a preferred direction i in reference configuration, and n is the total number of different preferred directions. The detailed definitions of structural tensors and the related theory of invariants have been reviewed by Zheng (1994)  .
Polyconvexity provides an alternative way to define invariant components. The polyconvex functions, , , and , describe deformations of distance, area, and volume elements since they map the corresponding elements from reference configuration to current configuration, respectively. Thus, the polyconvex functions play an important role in the definition of polyconvex invariant components. The invariant component defined in (12) happens to be polyconvex. Following Schröder and Neff (2003)  , another polyconvex invariant component in a preferred direction i, , can be defined as
The polyconvex invariant components, and , along with the invariant , have been selected as the arguments of an anisotropic stored energy functional by Itskov, Ehret, and Mavrilas (2006) 
Soft tissues can mainly be characterized as transversely isotropic and orthotropic hyperelastic materials. Transverse isotropy represents a material symmetry with respect to only one principal direction and the structural tensors can be expressed in the reference configuration as
where the principal material direction is denoted by the subscript 1. Orthotropy is characterized by symmetry with respect to three mutually orthogonal planes and the structural tensors read
Many influential discoveries in constitutive modeling of SBTs, including the structural tensors by Spencer (1971)  , the exponential strain energy function by Fung (1973)  , the anisotropic constitutive model by Holzapfel, Gasser, and Ogden (2000)  , the extension of polyconvexity to invariant components by Schröder and Neff (2003)  , and the selection of invariant components by Itskov, Ehret, and Mavrilas (2006)  , have been achieved. Existing stored energy functionals for modeling SBTs are usually assumed to be divided into three parts: a purely volumetric part, a purely isochoric part due to isotropic deformations, and a purely isochoric part due to anisotropic deformations. Anisotropic constitutive models are essentially composed of ICGs. The ICGs of existing models, however, are phenomenologically determined and invariant components are incompletely defined and inconsistently selected. To overcome such shortcomings, a physically relevant and mathematically covariant anisotropic CSE functional is needed for modeling and predicting finite deformations of SBTs.
The CSE functional, originally developed for modeling isotropic polymers by Zhao (2016, 2017)   , will be generalized to model anisotropic hyperelastic materials such as SBTs. The main objectives, therefore, are to rigorously formulate the anisotropic CSE functional with the ICGs generally grouped by the Lie group and particularly grouped by differential geometry, to newly define the second invariant component to capture anisotropic transverse effects and shear coupling effects, and to accurately fit anisotropic CSE models with experimental data for SBTs.
2. Anisotropic Constitutive Modeling
2.1. Anisotropic CSE Functional
An anisotropic CSE functional can be conceptually expressed as
where the range for subscript i, , is omitted for the rest of the equations for simplification and the number of different preferred directions, n = 3, is used for both orthotropic and transversely isotropic materials. The first invariant component can be expressed as
and the second invariant component is newly defined as
For convenient subsequent derivations, one can additionally define
The first-order derivatives of invariant components with respect to the right Cauchy-Green tensor are given below
The second-order derivatives of invariant components with respect to the right Cauchy-Green tensor have been obtained as
For the anisotropic CSE functional (17), along with the invariant components (18) through (20), the second Piola-Kirchhoff stress tensor for anisotropic hyperelastic materials, , reads
Substituting the first order derivatives (21), (22), and (23) into (27)2 yields
and the Kirchhoff stress tensor, , can be converted from the second Piola-Kirchhoff stress tensor by the following push-forward operation
For isothermal processes, the general CSE functional for anisotropic hyperelastic materials at finite deformations, , is postulated to be balanced with its stress work done as
Substituting (28) into (30)1 or (29) into (30)2, simplifying, and rearranging produces the CSE PDE
With Lie group methods, the characteristic system for the PDE is
and taking its three first-integrals, , , and , the general solution is obtained and written as
where and are two arbitrary functions and the general solution defines a group of anisotropic CSE functionals.
The general solution (33) has two arbitrary functions to be determined for practical applications. For translational and rotational deformations, the first arbitrary function, , is selected as
For different powers of ellipsoidal deformations, the inverse power function is required. Thus, the second arbitrary function, , is chosen as
Substituting (34) and (35) into (33) gives the particular solution of the anisotropic CSE functional
where the four coefficients for a preferred fiber direction i, , , , and , are unknown constitutive constants to be determined by experimental tests of anisotropic hyperelastic materials.
2.2. Anisotropic Elasticity Tensors
Elasticity tensors are crucial in studying mathematical properties of the constitutive relations. The anisotropic fourth-order elasticity tensor in material description is defined as
Substituting (27)2 into (37) and taking derivatives yields the fourth-order anisotropic elasticity tensor in material description
Substituting (21) through (26) into (38) yields
where and the twelve parameters are defined as
The fourth-order anisotropic elasticity tensor in spatial description, similar to the isotropic elasticity tensor  , is defined as
With , the anisotropic elasticity tensors can be further simplified.
2.3. CSE Models for Different Deformation Modes
Commonly used experimental tests are uniaxial tension, biaxial tension, and triaxial shear tests. Nominal stress and stretch are usually recorded in experimental tests. With the CSE functional (36), nominal stress as a function of stretch can be established by
The three derivatives of the CSE functional (36) for incompressible SBTs are given below
Substituting (48), , and derivatives of invariant components with respect to normal stretches into (47) yields the incompressible stress-free (ISF) condition in reference configuration
Thus, the CSE constitutive models in three deformation modes can be derived based on the equations (47) through (49).
2.3.1. Uniaxial Tension Mode
The deformation of uniaxial tension can be modeled as
where and denote the Cartesian coordinates of a typical particle in reference and current configurations, respectively. With ,
the tensors for the uniaxial tension mode in (50) are and . The first principal invariant is obtained as . In the longitudinal uniaxial tension for orthotropic materials, the structural tensor is . The two invariant components are and . Substituting (48), invariant components, and their derivatives into (47) yields the longitudinal nominal stress, , in uniaxial tension tests for orthotropic materials
and by the same token for , the transverse nominal stress is written as
For transversely isotropic materials, the structural tensors are defined in (15). In longitudinal tension tests with cylindrical specimens, , the tensile stretch is with isotropic contraction of . With the equation (47), the longitudinal nominal stress is given by
and in transverse tension tests, , the tensile stretch is with anisotropic contraction and the shortening of fiber is not necessarily identical to . The transverse nominal stress is
With the ISF condition (49), the CSE models (51), (52), (53), and (54) can be further simplified for subsequent curve fittings.
2.3.2. Biaxial Tension Mode
The deformation of biaxial tension can be generally modeled by
where or is the amount of shear. The tensors and for cases in (55) are given by
where can be determined by incompressible condition
For orthotropic materials with , , , and related derivatives, the longitudinal and transverse Cauchy normal stresses due to fiber 1 can be obtained as
and for orthotropic materials with , , , and related derivatives, similarly, the longitudinal and transverse Cauchy normal stresses due to fiber 2 turns out to be
Combining with longitudinal normal stress components and in (59) and (61) and summing up transverse normal stress components and in (60) and (62) yields the general CSE models for biaxial tension tests of orthotropic materials
Without shear deformations, , in-plane normal stretches become in-plane principal stretches and normal stresses become principal stresses. Ignoring both shear deformations and shear coupling effects produces the simplified equations. For the structure tensor , the simplified longitudinal Cauchy principal stress for equibiaxial tension tests is
and for the structure tensor , the simplified transverse Cauchy principal stress is
With the ISF condition (49), the CSE models (63), (64), (65), and (66) can be further simplified in subsequent curve fitting processes.
2.3.3. Triaxial Shear Mode
The deformation of triaxial shear can be modeled as
where the angle of shear is . The tensors , , and for cases in (67) are worked out as
and with the right Cauchy-Green tensor (68)2 the three principal invariants are obtained as and . For orthotropic materials, the structural tensor is defined in (16).
For the longitudinal shear, the structural tensor is . The two invariant components are . With all three zero derivatives in (47), the longitudinal shear stress component is
For the perpendicular shear, the structural tensor is . The two invariant components are and . The two derivatives are and . Substituting the derivatives (48) for i = 2, the two invariant components and , and their derivatives into the stress Equation (47) produces the perpendicular shear stress component
For the transverse shear, the structural tensor is . The two invariant components are , , and their derivatives are and . Substituting the derivatives (48) for i = 3, and , and their derivatives into the stress Equation (47) gives the transverse shear stress component
Combining shear stress components in longitudinal shear (70), perpendicular shear (71), and transverse shear (72) results in the CSE model for triaxial shear tests
3. Fitting CSE Models to Testing Data
A self-developed graphics digitizer with MATLAB has been used to read out pairs of stress values and corresponding stretches from stress-stretch curves cited in following subsections.
3.1. Modeling for Rabbit Abdominal Skins in Uniaxial Tension
Rabbit abdominal skins are usually treated as incompressible orthotropic hyperelastic materials and they have been tested in constrained uniaxial tension and biaxial tension. In the constrained uniaxial tension tests, the longitudinal tension is conducted under the fixed transverse stretch while the transverse tension is conducted under the fixed longitudinal stretch by Lanir and Fung (1974)  . The norminal stress-stretch experimental data in longitudinal and transverse directions has been calculated from the corresponding force-stretch curves and initial areas for sample 36 provided by Tong and Fung (1976)  . The uniaxial tension experimental data of rabbit abdominal skins in both longitudinal or 1 and transverse or 2 directions have been used to fit the uniaxial tension models (51) and (52) and the two sets of constitutive constants have been extracted by an iterative least square method. The comparison between the CSE model and the uniaxial tension test data of rabbit abdominal skins is shown in Figure 2.
Figure 2. Comparison between CSE Models and uniaxial tension tests for rabbit abdominal skins.
3.2. Modeling for Porcine Liver Tissues in Uniaxial Tension
Porcine liver tissues, composed of liver lobules and connective tissues, are transversely isotropic with the principal axis along the direction of the lobule. Uniaxial tension and compression tests of porcine liver tissues, using cylindrical specimens, have been conducted by Chui et al. (2007)  . Uniaxial tension experimental data averaging over five porcine liver tissue specimens in both longitudinal and transverse directions have been used to fit the uniaxial tension models (53) and (54). Two sets of constitutive constants have been solved by an iterative least square method. The comparison between the anisotropic CSE model and the uniaxial tension test data of porcine liver tissues is shown in Figure 3.
3.3. Modeling for Human Myocardial Tissues in Biaxial Tension
Heart walls consist of three distinct layers: the endocardium, the myocardium, and the epicardium. The ventricular myocardium is the main functional tissue of the heart wall. Biaxial loading is closer to the actual loading condition of passive ventricular myocardial tissues than uniaxial loading. Biaxial tension tests of human passive ventricular myocardial tissues have been conducted by Sommer et al. (2015)  . The Cauchy stress-stretch equibiaxial tension experimental data in the mean-fiber direction and the cross-fiber direction has been averaged among 26 specimens, respectively. The averaged equibiaxial tension experimental data have been used to fit the biaxial tension models (65) and (66) and the two
Figure 3. Comparison between CSE models and uniaxial tension tests for porcine liver tissues.
sets of constitutive constants for both directions have been obtained by an iterative least square method. The comparison between the anisotropic CSE models and the averaged equibiaxial tension test data of human passive ventricular myocardial tissues is shown in Figure 4.
3.4. Modeling for Human Myocardial Tissues in Triaxial Shear
Triaxial shear tests for porcine myocardial tissues demonstrate the orthotropic behavior of myocardial tissues by Dokos et al. (2002)  . The passive myocardial tissue can therefore be modeled as an orthotropic material with three orthogonal directions: fiber, sheet, and normal. The fiber refers to muscle fiber, the sheet covers a network of collagen fibers, and the normal direction is perpendicular to both muscular and collagen fiber directions  . The triaxial shear tests for human myocardial tissues have also been conducted by Sommer et al. (2015)  . Human myocardial tissues have also been modeled as orthotropic materials with three orthogonal directions: fiber, sheet, and normal. The data for triaxial perpendicular shear (TP-Shear) to fiber direction are averaged between those of fiber-sheet mode and fiber-normal mode, the data for TP-Shear to sheet direction are averaged between those of sheet-fiber mode and sheet-normal mode, and the data for TP-Shear to normal direction are averaged between those of normal-fiber mode and normal-sheet mode. The triaxial shear experimental data of ventricular myocardial tissues in fiber, sheet, and normal directions have been used to fit the CSE model for triaxial shear (73) and three sets of constitutive constants have been resolved by an iterative least square
Figure 4. Comparison between CSE models and biaxial tension tests for human myocardium.
method. The comparison between the CSE model and the triaxial shear test data of human passive myocardial tissues is shown in Figure 5.
3.5. Constitutive Constants
In anisotropic CSE constitutive models, four constitutive constants are needed to describe anisotropic finite deformations in a preferred direction. Based on the ISF condition (49), one of the four constitutive constants can be treated as a dependent constant. The first constant, , is treated as a dependent constant since the last three constants, , , and , are contained in the anisotropic elasticity tensors (39) and (46).
The ISF condition may be applied during extraction of constitutive constants for curve fittings of uniaxial tension and biaxial tension tests. There are two methods for curve fitting experimental data with the ISF condition. In the first method, three unknowns, , , and , are to be solved by the iterative least square method and is eliminated from a model by the ISF condition (49). In the second method, four unknowns, , , , and , are to be solved by the iterative least square method and the ISF condition (49) is simply subtracted from a model.
The constitutive constants for modeling rabbit abdominal skins (RAS) and porcine liver tissues (PLT) in uniaxial tension, human myocardial tissues (HMT) in biaxial tension and triaxial shear are collected in Table 1. Extraction of constitutive constants for uniaxial tensions in longitudinal (L) and transverse (T) directions, biaxial tensions in the mean-fiber (MF) and the cross-fiber (CF)
Figure 5. Comparison between CSE model and triaxial shear tests for human myocardium.
Table 1. Constitutive constants for different SBTs.
directions was conducted by the first method. Extraction of constitutive constants for triaxial shear tests in fiber (F), sheet (S), and normal (N) directions was conducted directly from the model since the ISF condition is automatically satisfied.
In the triaxial shear model (73), the extracted second constitutive constants, however, are the sums of and. For the three orthogonal directions of fiber, sheet, and normal, they are, , and, respectively. The data for triaxial perpendicular shear of fiber are averaged between those of fiber-sheet mode and fiber-normal mode. As an example for the fiber direction, we have and. Similarly, the two other equations can be established for the sheet and normal directions. Combining equations together for three orthogonal directions yields
Solving linear Equation (74) simultaneously gives
4.1. Anisotropic CSE Functional
In the anisotropic CSE functional, and instead of and are named and used albeit since the summation of invariant components or results in the corresponding invariant or
Moreover, the use of defined in (19) rather than in (13) captures anisotropic transverse deformations and avoids the unnecessary calculation of inverse right Cauchy-Green tensor.
In the anisotropic CSE functional (36), the first term, , represents the work done of normal stress and translational deformation. The second term, , describes the work done of shear stress and rotational deformation. The third term, , captures the work done of stress for finite deformations with different anisotropic degrees of ellipsoidal deformations.
With the fourth-order anisotropic elasticity tensor (38), the most general fourth-order elasticity tensor for isotropic hyperelastic materials can be recovered by the mappings of, , , …, and
Substituting the first and second order derivatives of invariants (6), (7), (8), and (9) into (77), simplifying, and rearranging produces
where the eight parameters are defined by
where the isotropic CSE functional is given by
The three constitutive constants, , and are generally used for modeling finite deformations of isotropic hyperelastic materials. The two constitutive constants, and, determine the isotropic elasticity tensor (78) with the eight parameters (79) through (82).
4.2. Comments on Different Experimental Tests
In uniaxial tension tests for SBTs, not only does the stress-stretch data in the tensile direction need to be captured but also at least one other orthogonal stretch must be measured for incompressible SBTs since needs to be evaluated for (19) in the CSE model, and transverse deformation effects on longitudinal deformations are indispensable in uniaxial tension tests. The transverse deformation effects in modeling SBTs has been elaborated on by Latorre, Romero, and Montáns (2016)  . With uniaxial tension measurements in the direction, for example, we also need to capture the corresponding transverse stretch or and use the incompressible condition to determine the other stretch or, respectively since and are generally not equivalent due to orthotropy. In constitutive modeling of the transverse uniaxial tension test of porcine liver tissues, the tensile stretch is and the shortenings of fibers and cross-fibers are and, respectively. Considering the anisotropic compression test results, the anisotropic contraction is simulated with and the constitutive constants have been obtained and listed in the PLT-T row in Table 1. If the isotropic contraction of were assumed, the constitutive constants would be and in which they are quite different from the PLT-T row of results in Table 1. Thus, it is essential to measure both the longitudinal and transverse stretches in uniaxial tension tests for incompressible SBTs.
In biaxial tension tests, it is very difficult to maintain both uniform force distribution and uniform normal deformations. In general, the stress-stretch curves in biaxial tension tests from cruciform specimens to square specimens are not accurate. For biaxial tension tests of square specimens, the stress as a function of stretch is generally overestimated. The overestimation, the correction factor, and the inverse finite element method regarding biaxial tension tests have been studied by Nolan and McGarry (2016)  . In the CSE models for biaxial tension tests, shear coupling effects do exist, making curve fittings harder. Without simplifications, at least four arguments, , , , and, have to be measured, causing material characterizations to be more complex.
In triaxial shear tests, according to different orthogonal fiber reinforcement orientations, shear deformations have been classified as longitudinal shear, perpendicular shear, and transverse shear by Destrade, Horgan, and Murphy (2015)  . The three orthogonal directions, fiber, sheet, and normal, for myocardial tissues generate six deformation modes: fiber-sheet mode, fibernormal mode, sheet-normal mode, sheet-fiber mode, normal-fiber mode, and normal-sheet mode. In the CSE model for triaxial shear tests, the longitudinal shear does not contribute any anisotropy, only the perpendicular shear elongates fiber reinforcements, and the perpendicular shear and transverse shear are coupled together. Without shear coupling effects, there would be three different experimental curves in the six deformation modes. The triaxial shear tests in six deformation modes for both porcine myocardial tissues by Dokos et al. (2002)  and human myocardial tissues by Sommer et al. (2015)  have produced more than three different curves, indicating shear coupling effects as well.
An anisotropic CSE functional, for isothermal processes, is postulated to be balanced with its stress work done, constructing a PDE. The anisotropic CSE PDE with respect to the three invariant components, , , and, is generally solved by Lie group methods. A three-term particular solution, which is essentially composed of ICGs, is particularly grouped by differential geometry to capture the three fundamental deformations. In a preferred direction i, the
term monitors translational deformations, the term captures
curvatures of rotational deformations, and the term describes curvatures of different powers of ellipsoidal deformations.
The anisotropic CSE constitutive models for uniaxial tension, biaxial tension, and triaxial shear tests have been derived. The constitutive constants have been solved by an iterative least square method for the uniaxial tension tests of rabbit abdominal skins and porcine liver tissues, the biaxial tension, and the triaxial shear tests of human ventricular myocardial tissues. With the newly defined second invariant component, the anisotropic CSE models capture the transverse effects in uniaxial tension deformations and the coupling effects between the perpendicular shear and transverse shear in triaxial shear deformations.
For anisotropic CSE models, the first constitutive constant, , can be treated as a dependent constant due to the ISF condition if necessary and the last three constitutive constants, , , and, are independent constants, defining the anisotropic elasticity tensor. For isotropic CSE models, the three constitutive constants, , , and, are needed for modeling finite deformations of isotropic hyperelastic materials. The two constitutive constants, and, are required for the isotropic elasticity tensor.
The author’s heartfelt thanks go to Jim S.J. Chen and Kenneth B. Margulies for their collaboration in earlier studies of rat hearts through ex-vivo tests of diastolic pressure as a function of balloon volume. He would also like to express his gratitude to Jianming and Jiesi Zhao for their discussion and assistance.