G Chagnon

M Rebouah

D Favier

Hyperelastic Energy Densities for Soft Biological Tissues: A Review

Many soft tissues are naturally made of a matrix and fibres that present some privileged directions. They are known to support large reversible deformations. The mechanical behaviour of these tissues highlights different phenomena as hysteresis, stress softening or relaxation. A hyperelastic constitutive equation is typically the basis of the model that describes the behaviour of the material. The hyperelastic constitutive equation can be isotropic or anisotropic, it is generally expressed by means of strain components or strain invariants. This paper proposes a review of these constitutive equations.

Introduction

Soft tissues are composed of several layers, each one of these layers has different compositions. It is considered that four typical tissues exist: epithelial tissue, connective tissue, muscular tissue and neuronal tissue [156]. For the mechanical studies on soft tissues the connective tissues are often considered as the most important from a mechanical point of view [START_REF] Epstein | Isolation and characterisation of cnbr peptides of human [α 1 (iii)] 3 collagen and tissue distribution [α 1 (i)] 2 α 2 and [α 1 (iii)[END_REF]156,177]. They are composed of cells and of extra cellular matrix. The extra cellular matrix is composed of ground substance and of three types of fibres: collagen, reticular and elastic fibres. Collagen fibres are often considered as more important than others, particularly because of their large size, and represent most of the mechanical behaviour. The reticular fibres, which are thin collagen fibres with different chemical properties, allow creating ramifications with the collagen fibres. Finally the elastic fibres mainly composed of elastin present a purely elastic behaviour and are also linked to the collagen fibres. The elastic properties of soft tissues are mainly due to these fibres. Soft tissues are often able to support large deformations.

The first mechanical study of soft tissues started in 1687 with Bernoulli experiments on gut. The first constitutive equation was proposed in 1690 by Leibniz, before Bernoulli and Riccati proposed other equations [START_REF] Bell | The Experimental Foundations of Solid Mechanics, Mechanics of Solids[END_REF]. Since these works, many experimental studies have been performed. As an illustration, some experimental data can be found, not exhaustively, in the literature for arteries [209,262], aortic valve tissues [162], veins [START_REF] Alastrué | Experimental study and constitutive modelling of the passive mechanical properties of the ovine infrarenal vena cava tissue[END_REF], vaginal tissues [196], anterior malleolar ligament [START_REF] Cheng | Mechanical properties of anterior malleolar ligament from experimental measurement and material modeling analysis[END_REF], muscles [START_REF] Gras | Hyper-elastic properties of the human sternocleidomastoideus muscle in tension[END_REF], human trachea [254], cornea [235], skin [START_REF] Groves | An anisotropic, hyperelastic model for skin: experimental measurements, finite element modelling and identification of parameters for human and murine skin[END_REF] or gallbladder walls [143]... Even if many soft tissues are studied, the largest database in the literature concerns arteries.

Soft tissues present a complex behaviour with many non-linear phenomena as explained by different authors [118,124] as the time dependency [START_REF] Bischoff | A rheological network model for the continuum anisotropic and viscoelastic behavior of soft tissue[END_REF]202] or the stress softening phenomenon [154,198], i.e., their mechanical behaviour mainly depends on time and on the maximum deformation previously endured. Most of soft tissues dissipate energy when loading, nevertheless, the elastic behaviour generally dominates their behaviour and it represents the asymptotic behaviour when the dissipation diminishes to zero. In this way, in a first approach, most of the soft tissues are described in the context of hyperelasticity [START_REF] Harb | A new parameter identification method of soft biological tissue combining genetic algorithm with analytical optimization[END_REF]149,246]. To take into account the fibrous structure of the soft tissues, anisotropic formalism is introduced. The diversity among the mechanical characteristics of soft tissues has motivated a great number of constitutive formulations for the different tissue types. For example, the reader is referred to [224], wherein the author treats the history of biaxial techniques for soft planar tissues and the associated constitutive equations.

Anisotropic hyperelasticity can be modeled by using the components of the strain tensor or by the use of strain invariants. The two formulations permit the development of different families of anisotropic strain energy densities. Soft tissues are numerous and present different tissue architectures that lead to various anisotropy degrees, i.e., difference of mechanical behaviour in each direction, and different maximum admissible deformation. In this way, many constitutive equations are proposed to describe the tissues.

The aim of this paper is to propose a review of most of the hyperelastic strain energy densities commonly used to describe soft tissues. In a first part, the different formalisms that can be used are recalled. In a second part, the isotropic modelling is described. In a third part, the anisotropic modelling is presented. The deformation tensor component approach based on Fung's formulation is briefly presented, and invariant approaches are detailed. In a fourth part, the statistical approaches, considering the evolution of the collagen network, are described. Last, a discussion about the models closes the paper.

Mechanical Formulation

Description of the Deformation

Deformations of a material are classically characterised by right and left Cauchy-Green tensors defined as C = F T F and B = FF T , where F is the deformation gradient. In the polar decomposition of F, the principle components of the right or left stretch tensors are called the stretches and are denoted as λ i with i = 1..3. The Green-Lagrange tensor is defined as E = (C -I)/2, where I is the identity tensor, and its components are denoted as E ij with i, j = 1...3. Nevertheless, some prefer to use the logarithmic strains e i = ln(λ i ), instead of a strain tensor, generalised strains as e i = 1 n (λ n i -1) [185], or others measures as, for example,

e i = λ i λ 2 j λ 2 k
with j = i and k = i [START_REF] Gilchrist | Generalisations of the strain-energy function of linear elasticity to model biological soft tissue[END_REF]; all these measures are written in their principal basis.

Instead of using directly the strain tensors, strain invariants are often preferred as they have the same values whatever the basis is. From an isotropic point of view, three principal strain invariants I 1 , I 2 and I 3 are defined by

I 1 = tr(C), (1)
I 2 = 1 2 tr(C) 2 -tr C 2 , ( 2
)
I 3 = det(C), (3)
where "tr" is the trace operator, and "det" the determinant operator. Characteristic directions corresponding to the fibre orientations must be defined. For one material, one or many material directions (the number of directions is noted q) can be defined according to the architecture of the considered tissue. In the undeformed state, the ith direction is noted N i in the initial configuration. The norm of the vector N i is unit. Due to material deformation, the fibre orientations are evolving in the deformed state. The current orientation is defined by n (i) = FN (i) .

(4)

Note that n (i) is not a unit vector. Two orientation tensors can be defined, one in the undeformed and another in the deformed state:

A (i) = N (i) ⊗ N (i) , a (i) = n (i) ⊗ n (i) . ( 5
)
The introduction of such directions lead to the definition of new invariants related to each direction. The invariant formulation of anisotropic constitutive equations is based on the concept of structural tensors [START_REF] Boehler | Applications of Tensor Functions in Solid Mechanics[END_REF][START_REF] Boehler | A simple derivation of representations for non-polynomial constitutive equations in some cases of anisotropy[END_REF]238,239,241,243]. 1 The invariant I 4 and I 5 can be defined for one direction i as

I (i) 4 = tr CA (i) = N (i)
• CN (i) , and

I (i) 5 = tr C 2 A (i) = N (i) • C 2 N (i) . ( 6
)
In practice, some prefer to use the cofactor tensor of F, i.e., Cof(F), [120] and to define J (i) 5 = tr(Cof(C)A (i) ), in order to easily ensure the polyconvexity of the strain energy (see Sect. 2.3). In the literature, in the case of two fibre directions (1) and (2), a notation I 4 and I 6 is often used for soft tissues [108] instead of I (1) 4 and I (2) 4 (or I 5 and I 7 instead of I (1) 5 and I (2) 5 ). In this paper, it is preferred to keep only I (i) 4 notation and to generalise the notation to n directions. These invariants depend only on one direction but it is possible to take into account the interaction between different directions, by introducing a coupling between directions i and j by means of two other invariants:

I (i,j ) 8
= N (i) • N (j ) N (i) • CN (j ) , and

I (i,j ) 9 = N (i) • N (j ) 2 . ( 7
)
I (i,j ) 9
is constant during deformation, thus it is not adapted to describe the deformation of the material but it represents the value of I (i,j ) 8 for zero deformation. Let us denote I k as the invariants family (I 1 , I 2 , I 3 , I (i) 4 , I (i) 5 , I (i,j ) 8

, I (i,j ) 9

) and J k as the invariants family 1 Details about the link between structural tensors and a method to link a fictitious isotropic configuration to render an anisotropic, undeformed reference configuration via an appropriate linear tangent map is given in [163].

(I 1 , I 2 , I 3 , I (i) 4 , J (i) 5 ). When only one direction is considered, the superscript (i) is omitted in the remainder of this paper.

The I k invariants are the mostly used invariants in the literature, although other invariants have been proposed. Some authors [START_REF] Ciarletta | Stiffening by fiber reinforcement in soft materials: a hyperelastic theory at large strains and its application[END_REF] propose to use invariants that are zero at zero deformation. In this way, they introduce the tensor G = H T H, with H = 1 2 (F -F T ). This motivates the definition of a new class of invariants

I k : ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ I 1 = tr(G), I 2 = tr G 2 , I 4 = tr GA (i) , I 5 = tr G 2 A (i) . ( 8
)
Ericksen and Rivlin [START_REF] Ericksen | Large elastic deformations of homogeneous anisotropic materials[END_REF] proposed another formulation, adapted to transversely isotropic materials only, characterised by a vector N (i.e., only one direction i). This direction often corresponds to a fibre reinforced direction. Their work was further used by different authors [START_REF] Agoras | A general hyperelastic model for incompressible fiber-reinforced elastomers[END_REF][START_REF] Criscione | Physically based strain invariants set for materials exhibiting transversely isotropic behavior[END_REF][START_REF] Debotton | Neo-Hookean fiber-reinforced composites in finite elasticity[END_REF][START_REF] Debotton | Mechanics of composites with two families of finitely extensible fibers undergoing large deformations[END_REF] who proposed to define other invariants (λ p , λ n , γ n , γ p , ψ γ ), denoted as Cr k . They can be expressed as a function of the I k invariants:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λ 2 p = I 3 I 4 , λ 2 n = I 4 , γ 2 n = I 5 I 4 -I 4 , γ 2 p = I 1 - I 5 I 4 -2 I 3 I 4 , tan 2ψ γ = 2λ p H + / -γ p γ 4 n γ 2 p (4λ 2 p + γ 2 p ) -H 2 λ p H + / -2λ p γ 4 n γ 2 p (4λ 2 p + γ 2 p ) -H 2 , ( 9
)
with

H = (2λ 2 n + γ 2 n )(2λ 2 p + γ 2 p ) + 2λ 4 p -2I 2 .
The advantage is that these invariants have a physical meaning. λ n is the measure of stretch along N, λ p is a measure of the in-plane transverse dilatation, γ n is a measure of the amount of out-of-plane shear, γ p is the amount of shear in the transverse plane, and ψ γ is a measure of the coupling among the other invariants. Criscione et al. [START_REF] Criscione | Physically based strain invariants set for materials exhibiting transversely isotropic behavior[END_REF] criticised these invariants for not being zero for zero deformation, as is the corresponding strain tensors. They proposed to use the β k invariants:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ β 1 = ln I 3 2 , β 2 = 3 ln I 4 -ln I 3 4 , β 3 = ln I 1 I 4 -I 5 2 √ I 3 I 4 + I 1 I 4 -I 5 2 √ I 3 I 4 2 -1 , β 4 = I 5 I 2 4
-1,

β 5 = I 1 I 4 I 5 + I 1 I 3 4 + 2I 3 I 4 -I 2 5 -2I 2 I 2 4 -I 5 I 2 4 (I 5 -I 2 4 ) I 2 1 I 2 4 + I 2 5 -2I 1 I 4 I 5 -4I 3 I 4 . ( 10
)
These invariants also have a physical meaning. β 1 is the logarithmic volume strain, β 2 specifies a fibre strain of distortion, β 3 specifies the magnitude of cross-fibre, i.e., pure shear strain, β 4 specifies the magnitude of along fibre strain, i.e., simple shear strain and β 5 specifies the orientation of the along fibre shear strain relative to the cross-fibre shear strain. These last two families of invariants were developed for a one fibre direction material; it can easily be generalised to q directions but this has not yet been used yet in the literature. All these invariants are useful. In practice, the I k are the most used and the other invariants are not often used for calculations in finite element software. But, as they can be written by means of the I k invariants, all the expressions can be deduced from these invariants. As a consequence in this work, the theoretical development is only presented for the I k formulation.

Strain-Stress Relationships

Living tissues are often considered as incompressible. To use constitutive equations in finite element codes, a volumetric/isochoric decomposition is used. All the equations are written using the pure incompressibility hypothesis in order to avoid any non-physical response of these equations [100], but some details about the consequences of the volumetric-isochoric choice split is detailed in [227]. Nevertheless, they can be written in a quasi-incompressible framework by means of the incompressible invariants Īk :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ī1 = I 1 I -1/3 3 , Ī2 = I 2 I -2/3 3 , Ī4 (i) = I (i) 4 I -1/3 3 , Ī5 (i) = I (i) 5 I -2/3 3 , Ī8 (i,j ) = I (i,j ) 8 I -1/3 3 . ( 11
)
This formulation is particularly useful for finite element implementation. All the equations for the elasticity tensor can be seen in different papers [START_REF] Bose | Computational aspects of a pseudo-elastic constitutive model for muscle properties in a soft-bodied arthropod[END_REF]133,145,152,195,271]. In this case, a penalty function depending on I 3 is used to ensure incompressibility. One can refer to [START_REF] Doll | On the development of volumetric strain energy functions[END_REF] for a comparison of the different functions classically used. The choice of the penalty parameter to ensure incompressibility [253] is a critical issue. In this paper, all the constitutive equations are written in the purely incompressible framework, but all the models can be established in the quasi-incompressible framework as well.

The second Piola-Kirchhoff stress tensor can be directly calculated by derivation of the strain energy function W (I 1 , I 2 , I (i) 4 , I (i) 5 , I (i,j ) 8

, I (i,j ) 9

), with i, j = 1..q:

S = 2 (W ,1 + I 1 W ,2 )I -W ,2 C + q i W (i) ,4 N (i) ⊗ N (i) + q i W (i) ,5 N (i) ⊗ CN (i) + N (i) C ⊗ N (i) + i =j W (i,j ) ,8 N (i) • N (j ) N (i) ⊗ N (j ) + N (j ) ⊗ N (i) + pC -1 (12)
where W ,k = ∂W ∂I k , and p is the hydrostatic pressure. The Eulerian stress tensor, i.e., the Cauchy stress tensor, is directly obtained by the push-forward operation. To ensure that the stress is identically zero in the undeformed configuration, it is required that: ∀i W (i) ,4 + 2W (i) ,5 = 0, [START_REF] Baek | Theory of small on large: potential utility in computations of fluid-solid interactions in arteries[END_REF] for zero deformation [174]. The direct expressions that permit calculation of the stress with the other invariants basis can be found in [START_REF] Criscione | Physically based strain invariants set for materials exhibiting transversely isotropic behavior[END_REF][START_REF] Criscione | Constitutive framework optimized for myocardium and other high-strain, laminar materials with one fiber family[END_REF].

Stability

The strong ellipticity condition is a mathematical restriction on the constitutive functions.

For three-dimensional problems [267], the strong ellipticity was characterised for compressible isotropic materials in [236], and for incompressible ones in [273]. In this context, the strong ellipticity was largely studied in the case of transverse isotropy for in plane strains in [165-168, 170, 240]. The generic condition to verify for the strain energy in the absence of body forces [150,167,168] can be written as:

1 J F pr F qs ∂ 2 W ∂F ir F js n p n q m i m j > 0 with m = 0 and n = 0, (14)
where m and n are two non-zero vectors. Nevertheless, this condition is always difficult to verify. Thus, some have proposed another way to tackle the strong ellipticity condition.

It is known that polyconvexity implies ellipticity [173,228,232]. As a consequence, the polyconvexity in the sense of Ball [START_REF] Ball | Convexity conditions and existence theorems in non-linear elasticity[END_REF][START_REF] Ball | Constitutive equalities and existence theorems in elasticity[END_REF] is used, even if it is more restrictive than strong ellipticity. Of course, some strain energies can be elliptic but not polyconvex. It is important to note that polyconvexity does not conflict with the possible non-uniqueness of equilibrium solutions, as it guarantees only the existence of at least one minimizing deformation. Hence, polyconvexity provides an excellent starting point to formulate strain energy functions that guarantees both ellipticity and existence of a global minimizer.

Polyconvexity has been studied within the framework of isotropy [START_REF] Bilgili | Restricting the hyperelastic models for elastomers based on some thermodynamical, mechanical and empirical criteria[END_REF]244], and the conditions to verify it are well known for every classical isotropic model from the literature (see for example [START_REF] Hartmann | Parameter estimation of hyperelasticity relations of generalized polynomial-type with constraint conditions[END_REF][START_REF] Hartmann | Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility[END_REF]180,204]). Many authors have extended their study to anisotropic materials [START_REF] Ehret | A polyconvex hyperelastic model for fiber-reinforced materials in application to soft tissues[END_REF]121,171,206,245,257]. Some have studied the polyconvexity of existing constitutive equations [START_REF] Doyle | Adaptation of a rabbit myocardium material model for use in a canine left ventricle simulation study[END_REF]104,106,186,267], whereas others have attempted to directly develop polyconvex constitutive equations. Some Conditions. In case of existing constitutive equations, Walton and Wilber [267] summarised conditions to ensure polyconvexity. For a strain energy depending on I 1 , I 2 and I 4 , W (I 1 , I 2 , I (i) 4 ), the conditions are:

W ,k > 0 for k = 1, 2, 4 and (15)
[W ,kl ] is definite positive. ( 16)

If the strain energy also depends on

I (i,j ) 8
, the following condition should be added:

∂W ∂I (i,j ) 8 ≤ ∂W ∂I 1 . ( 17
)
The use of the fifth invariant I (i) 5 introduces the need to change the other invariants, as I 5 is not a polyconvex function (when used alone). Walton and Wilber [267]

used I * k : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ I * 1 = 1 2 I 1 , I * 2 = 1 2 I 2 1 -I 2 , I (i) * 4 = I (i) 4 , I (i) * 5 = I (i) 5 . ( 18
)
Here, the condition to verify for

I * k is: ⎧ ⎪ ⎨ ⎪ ⎩ W ,k > 0 for k = 1, 2, 4, 5, W ,1 + κW ,4 ≥ 0 for some κ > 4, [W ,kl ] is definite positive. ( 19
)
As it will be described in next paragraph, many strain energies can be decomposed as W = W iso (I 1 ) + W aniso (I 4 ). In this case, some sufficient conditions, but not necessary for polyconvexity, have been given in [106] for the anisotropic part:

∂W aniso ∂I 4 ≥ 0 and ( 20)

∂W aniso ∂I 4 + 2I 4 ∂ 2 W aniso ∂I 2 4 ≥ 0. (21)
These two restrictive conditions mean that the considered directions cannot generate negative forces when submitted to compression whereas the strong ellipticity can also be verified in compression. This is an illustration of the constraints generated by the polyconvexity compared to strong ellipticity.

Development of Specific Constitutive Equations. Some authors have created elementary strain energies that satisfy polyconvexity. First, Schroder and Neff [228] worked on equations depending on I 1 and I 4 , and they proved that some functions are polyconvex:

W 1 = β 1 I 4 , W 2 = β 2 I 2 4 , W 3 = β 3 I 4 I 1/3 3
, and

W 4 = β 4 I 2 4 I 1/3 3 , ( 22
)
where β i are material parameters. Nevertheless, as I 5 is not a polyconvex function, some have proposed [228,229] the construction of new combinations of invariants in the case of one reinforced direction that are polyconvex; these invariants are denoted as

K i : ⎧ ⎪ ⎨ ⎪ ⎩ K 1 = I 5 -I 1 I 4 + I 2 tr(A) 1/2 , K 2 = I 1 -I 4 , K 3 = I 1 I 4 -I 5 . ( 23
)
These invariants permitted the development of a list of elementary polyconvex energies [START_REF] Ebbing | Approximation of anisotropic elasticity tensors at the reference state with polyconvex energies[END_REF]231]. The different strain energies are listed in Table 1. Since a combination of polyconvex energy densities is also polyconvex, it is possible to develop many constitutive equations that can be adapted to different soft tissues.

Table 1 Elementary polyconvex functions [START_REF] Ebbing | Approximation of anisotropic elasticity tensors at the reference state with polyconvex energies[END_REF]231], where β i with i = 5...23 are material parameters Elementary polyconvex functions

W 5 = β 5 K 1 W 6 = β 6 K 2 1 W 7 = β 7 K 3 1 W 8 = β 8 K 1 I 1/3 3 W 9 = β 9 K 2 1 I 2/3 3 W 10 = β 10 K 2 W 11 = β 11 K 2 2 W 12 = β 12 K 2 I 1/3 3 W 13 = β 13 K 2 2 I 2/3 3 W 14 = β 14 K 3 W 15 = β 15 K 2 3 W 16 = β 16 K 3 I 1/3 3 W 17 = β 17 K 2 3 I 2/3 3 W 18 = β 18 I 2 1 + I 4 I 1 W 19 = β 19 2I 2 2 + I 2 I 5 -I 1 I 2 I 4 W 20 = β 20 3I 2 1 -I 4 I 1 W 21 = β 21 2I 2 2 + I 1 I 2 I 4 -I 2 I 5 W 22 = β 22 (3I 1 -2I 4 ) W 23 = β 23 (I 2 -2I 5 + 2I 1 I 4 )

Isotropic Hyperelastic Constitutive Equations

From a macroscopic point of view, soft tissues are an assembly of cells and fibres. According to the quantity and the orientation of the fibres, the behaviour of soft tissues can be supposed isotropic or not. According to the application, anisotropic behaviour can be neglected, and isotropic modelling can be efficient. In this way, many authors decide to use an isotropic approach to model soft tissues, as for example liver [149] kidney [113], bladder and rectum [START_REF] Boubaker | Finite element simulation of interactions between pelvic organs: predictive model of the prostate motion in the context of radiotherapy[END_REF], pelvic floor [193], breast [START_REF] Azar | A deformable finite element model of the breast for predicting mechanical deformations under external perturbations[END_REF]226], cartilage [144], meniscus [START_REF] Abraham | Hyperelastic properties of human meniscal attachments[END_REF], ligaments [START_REF] Garcia | A nonlinear biphasic viscohyperelastic model for articular cartilage[END_REF], eardrum [START_REF] Cheng | Viscoelastic properties of human tympanic membrane[END_REF], arteries [192], brain [127], lungs [234], uterus [START_REF] Harrison | Towards a novel tensile elastometer for soft tissue[END_REF] or skin [142]... Many models that are used to describe an isotropic approach come from rubber like materials studies. Some literature reviews have been proposed [START_REF] Boyce | Constitutive models of rubber elasticity: a review[END_REF]258]. Constitutive equations for rubber like materials were created to represent a strain hardening for deformations of about hundreds of percent whereas soft tissues often strain harden after some tens of percent. Thus, the functions for rubber like materials may not necessarily apply. Other, more suitable constitutive equations have been developed especially for soft tissues. The main models are listed in Table 2. The main feature for the constitutive equations is the presence of an important change of slope in the strain-stress curve for moderate deformations. This explains why most of the equations include an exponential form which allows the description of strong slope changes. Nevertheless, all constitutive equations stay equivalent to the neo-Hookean model [255,256] for small strains. Moreover, most of the constitutive equations are very similar for the I 1 part as it is the exponential form that dominates in the equations. While most of the constitutive equations are only expressed with the first invariant, the second invariant can be employed to capture the different states of loading [112]. There exists some limitations to use only the first invariant [110,270]. Nevertheless, the choice of using I 1 , or (I 1 , I 2 ) mainly depends on the available experimental data. When experiments are limited to one loading case, it can be difficult to correctly fit a constitutive equation expressed by means of the two invariants.

Anisotropic Hyperelastic Constitutive Equations

Different approaches have been used to describe the anisotropy of soft tissues. The first one is based on Green-Lagrange components and the second one is based on strain invariants.

W = c 1 (I 1 -3) + c 2 (I 1 -3) 2 Knowles [131, 274] (*) W = c 1 2c 2 1 + c 2 c 3 (I 1 -3) c 3 -1 Exponential model Demiray [57] (**) W = c 1 c 2 exp c 2 2 (I 1 -3) -1 Demiray et al. [58] W = c 1 c 2 exp c 2 2 (I 1 -3) 2 -1 Holmes and Wow [103] W = c 0 exp c 1 (I 1 -3) + exp c 2 (I 2 -3) -c 0 Arnoux et al. [7, 8] W = c 1 exp c 2 (I 1 -3) - c 1 c 2 2 (I 2 -3) Singh et al. [237] W = c 1 2c 2 exp c 2 (I 1 -3) -1 + c 3 2 (I 2 -3) 2 Volokh and Vorp [266] W = c 1 -c 1 exp - c 2 c 1 (I 1 -3) - c 3 c 1 (I 1 -3) 2 Tang et al. [251] W = c 1 (I 1 -3) + c 2 (I 2 -3) + c 3 exp c 4 (I 1 -3) -1 Van Dam et al. [261] W = c 1 - 1-c 2 c 2 3 (c 3 x + 1) exp(-c 3 x) -1 + 1 2 c 2 x 2 with x = √ c 4 I 1 + (1 -c 4 )I 2 -3

Use of Green-Lagrange Tensor Components

The first model using the components of the Green-Lagrange strain tensor were developed in [118]. It consists in proposing strain energy densities that are summarily decomposed into contributions of each component with different weights; a review of these models is proposed in [116]. The first generic form was proposed by Tong and Fung [252]:

W = c 2 exp b 1 E 2 11 + b 2 E 2 22 + b 3 E 2 12 + E 2 21 + 2b 4 E 12 E 21 + b 5 E 3 11 + b 6 E 3 22 + b 7 E 2 11 E 22 + b 8 E 11 E 2 22 -1 , ( 24
)
where c and b i , i = 1...8 are material parameters. Three years later, Fung [START_REF] Fung | Pseudoelasticity of arteries and the choice of its mathematical expression[END_REF] developed a generic form in two dimensions, the model was next generalised to three dimensions [START_REF] Chuong | Three-dimensional stress distribution in arteries[END_REF].

Later, shear strains were introduced [128], and finally a global formulation was proposed [116]:

W = c exp(A ij kl E ij E kl ) -1 , ( 25
)
where c and A ij kl are material parameters. Different constitutive equations were then developed and written in cylindrical coordinates (r, θ , z) often used for arteries [138]. Moreover, the strain energy function can be naturally uncoupled into a dilatational and a distortional part [START_REF] Ateshian | A frame-invariant formulation of fung elasticity[END_REF], to facilitate the computational implementation of incompressibility. In the same way, as in non-Gaussian theory [137], it is possible to take into account the limiting extensibility of the fibres [175]. This exposes the possibility of a constitutive equation that presents an asymptote even if constitutive equations that include an exponential or an asymptotic form can be very close [START_REF] Chagnon | A comparison of the physical model of Arruda-Boyce with the empirical Hart-Smith model and the Gent model[END_REF]. The proposed models are listed in Table 3. The main difficulty of these constitutive equations is that they have a large number of material parameters.

Q = A ij kl E ij E kl Fung et al. [77] Q = b 1 E 2 θθ + b 2 E 2 zz + 2b 4 E θθ E zz Chuong and Fung [49] Q = b 1 E 2 θθ + b 2 E 2 zz + b 3 E 2 rr + 2b 4 E θθ E zz + 2b 5 E rr E zz + 2b 6 E θθ E rr Humphrey [116] Q = b 1 E 2 θθ + b 2 E 2 zz + b 3 E 2 rr + 2b 4 E θθ E zz + 2b 5 E rr E zz + 2b 6 E θθ E rr + b 7 E 2 θz + b 8 E 2 rz + b 9 E 2 rθ Costa et al. [51] Q = b 1 E 2 ff + b 2 E 2 ss + b 3 E 2 nn + 2b 4 1 2 (E f n + E nf ) 2 + 2b 5 1 2 (E sn + E ns ) 2 + 2b 6 1 2 (E f s + E sf ) 2 Rajagopal et al. [213] Q = b 1 E 2 θθ + b 2 E 2 zz + b 3 E 2 rr + 2b 4 E θθ E zz + 2b 5 E rr E zz + 2b 6 E θθ E rr + b 7 E 2 rr + E 2 θθ + b 8 E 2 θθ + E 2 zz + b 9 E 2 rr + E 2 zz
Other exponential functions Choi and Vito [START_REF] Choi | Two-dimensional stress-strain relationship for canine pericardium[END_REF]

W = b 0 exp b 1 E 2 11 + exp b 2 E 2 22 + exp(2b 3 E 11 E 22 ) -3 Kasyanov and Rachev [128] W = b 1 exp b 2 E 2 zz + b 3 E zz E θθ + b 4 E 2 θθ + b 5 E 2 zz E θθ + b 6 E zz E 2 θθ -1 + b 7 E θθ exp(b 8 E θθ ) + b 9 E zz + b 10 E 2 θz Other models Vaishnav et al. [259] W = b 1 E 2 θθ + b 2 E θθ E zz + b 3 E 2 zz + b 4 E 3 θθ + b 5 E 2 θθ E zz + b 6 E θθ E
Humphrey [117] W = b 1 E 2 rr + b 2 E 2 θθ + b 3 E 2 zz + 2b 4 E rr E θθ + 2b 5 E θθ E zz + 2b 6 E rr E zz + b 7 E 2 rθ + E 2 θr + b 8 E 2 zθ + E 2 θz + b 9 E 2 zr + E 2 rz Takamizawa and Hayashi [249] W = -c ln 1 -1 2 b 1 E 2 θθ + 1 2 b 2 E 2 zz + b 3 E θθ E zz + b 4 E θθ E zz + b 5 E θθ E rr + b 6 E rr E zz

Use of Strain Invariants

Strain energy densities depend on isotropic and anisotropic strain invariants. The use of I 4 and I 5 is necessary to recover linear theory [174]. Different cases exist. In a first case, the strain energy can be split as a sum into different parts as an isotropic and anisotropic contribution:

W = W iso (I 1 , I 2 ) + i W aniso I (i) 4 , I (i) 5 , ( 26
)
or some coupling can be realised between the isotropic and anisotropic parts as W aniso (I 1 , I 2 , I (i) 4 , I (i) 5 ). But very few models present a non-additive decomposition between two directions i and j , i.e., between

I (i) 4 , I (i) 5 , I (j )
4 and I (j )

5 . When W iso is used, it is often represented by a classical energy function. We discuss W aniso in the next paragraph. The use of only I 4 or I 5 , instead of the both of these invariants is questionable as it leads to the same shear modulus in the direction of and in the direction orthogonal to the reinforced direction [174].

Different model forms can be distinguished such as the polynomial, the power, the exponential and other constitutive equations not of these types.

Polynomial Development

The most known model for isotropic hyperelasticity is Rivlin's series [217] that describes a general form of constitutive equations depending on the first and second invariants. The generalisation of this model to an anisotropic formulation has been proposed in different ways. One consists in introducing the anisotropic invariants in the series. First a simple I 4 series [123] was proposed:

W aniso = n k=2 c i (I 4 -1) k , ( 27
)
where c i are material parameters. A linear term cannot be used, i.e., k = 1 in the previous equation, as it does not ensure zero stress for zero deformation. The term k = 2 corresponds to the standard reinforcing model [START_REF] Destrade | Surface instability of sheared soft tissues[END_REF]182,220,257], not initially proposed for soft tissues.

The complete generalisation of the Rivlin series was proposed in [222]:

W = klmn c klmn (I 1 -3) k (I 2 -3) l (I 4 -1) m (I 5 -1) n , ( 28
)
where c klmn are material parameters. A modified formulation was proposed in [111] to be more convenient for numerical use:

W = klmn c klmn (I 1 -3) k (I 2 -3) -3(I 1 -3) l (I 4 -1) m (I 5 -2I 4 + 1) n . ( 29
)
Instead of using I 4 , one may use √ I 4 which represents the elongation in the considered direction. This leads to a new series development [115]: W aniso = c 2 (I 4 -1) 2 + c 4 (I 4 -1) 4 Basciano and Kleinstreuer [START_REF] Basciano | Invariant-based anisotropic constitutive models of the healthy and aneurysmal abdominal aortic wall[END_REF] W aniso = c 2 (I 4 -1) 2 + c 3 (I 4 -1) 3 + c 4 (I 4 -1) 4 + c 5 (I 4 -1) 5 + c 6 (I 4 -1) 6 Basciano and Kleinstreuer [START_REF] Basciano | Invariant-based anisotropic constitutive models of the healthy and aneurysmal abdominal aortic wall[END_REF] W aniso = c 6 (I 4 -1) 6

W = kl c kl (I 1 -3) k ( I 4 -1) l , ( 30
)

Lin and Yin

[148] W = c 1 (I 1 -3)(I 4 -1) + c 2 (I 1 -3) 2 + c 3 (I 4 -1) 2 + c 4 (I 1 -3) + c 5 (I 4 -1) √ I 4 forms Alastrue et al. [4, 6] W aniso = c 2 ( √ I 4 -1) 2 Humphrey [115] W = c 1 ( √ I 4 -1) 2 + c 2 ( √ I 4 -1) 3 + c 3 (I 1 -3) + c 4 (I 1 -3)( √ I 4 -1) + c 5 (I 1 -3) 2 I 4 , I 5 forms

Park and Youn [194]

W aniso = c 3 (I 4 -1) + c 5 (I 5 -1) Bonet and Burton [START_REF] Bonet | A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations[END_REF] W Ogden [111,169,170] W aniso = c 2 (I 5 -1) 2

= c 1 + c 2 (I 1 -3) + c 3 (I 4 -1) (I 4 -1) - c 1 2 (I 5 -1) Bonet and Burton [31] W aniso = c 1 + c 3 (I 4 -1) (I 4 -1) - c 1 2 (I 5 -1) Merodio and

Hollingsworth and Wagner [102]

W aniso = c 2 (I 5 -

I 2 4 ) Murphy [174] W = c 1 (I 1 -3) + c 2 (2I 4 -I 5 -1) + c 3 (I 5 -1) 2 Murphy [174] W = c 1 (I 1 -3) + c 2 (2I 4 -I 5 -1) + c 3 (I 4 -1)(I 5 -1) Murphy [174] W = c 1 (I 1 -3) + c 2 (2I 4 -I 5 -1) + c 3 (I 4 -1) 2
where c kl are material parameters. It is worth noting that the use of √ I 4 includes, in the quadratic formulation [START_REF] Alastrué | On the use of the bingham statistical distribution in microsphere-based constitutive models for arterial tissue[END_REF][START_REF] Brown | A simple transversely isotropic hyperelastic constitutive model suitable for finite element analysis of fiber reinforced elastomers[END_REF], a model that represents the behaviour of a linear spring.

As other invariants were proposed, a series development based on β k invariants also has been considered [START_REF] Criscione | Physically based strain invariants set for materials exhibiting transversely isotropic behavior[END_REF]:

W = klm G klm β k 3 β l 4 β m 5 , ( 31
)
where G klm are material parameters.

As for rubber like materials with the Rivlin's series, the whole series is not used and a good truncation of the strain energy is essential. According to the considered material and to the loading states, different developments have been given in the literature. A list of resulting equations is included in Table 4. It is also important to note that the I 4 invariant is often used whereas the I 5 invariant is most often disregarded.

Power Development

Ogden 's [184] isotropic constitutive equation has proved its efficiency to describe complex behaviour. It is based on elongations and a power law development. For a material with a Schroder et al. [START_REF] Ciarletta | Stiffening by fiber reinforcement in soft materials: a hyperelastic theory at large strains and its application[END_REF]228,230]

W aniso = k 1 I k 2 4
Balzani et al. [START_REF] Balzani | A polyconvex framework for soft biological tissues. Adjustement to experimental data[END_REF][START_REF] Balzani | Simulation of discontinuous damage incorporating residual stresses in circumferentially overstretched atherosclerotic arteries[END_REF]

W = k 1 (I 1 I 4 -I 5 -2) k 2 Schroder et al. [230] W = k 1 (I 5 -I 1 I 4 + I 2 ) + k 2 I k 3 4 + k 4 (I 1 I 4 -I 5 ) + k 5 I k 6 4 O'Connell et al. [183] W = k 6 I 4 (I 5 -I 1 I 4 + I 2 ) -1 2
single fibre direction, there is the following generic form [186]:

W aniso = 2μ 1 β 2 I β/2 4 + 2I -β/4 4 -3 , ( 32
)
where μ 1 and β are material parameters. A generalised form was proposed by not imposing the same parameters for the two terms [264]

W aniso = r α r I βr 4 -1 + γ r I -δr 4 -1 , (33)
where α r , β r , γ r , δ r are material parameters. The same type of formulation is also proposed using the other invariants. Two general equations are of the form [122]:

W = klmn c klmn (I 1 -3) a k (I 2 -3) b l (I 4 -1
) cm (I 5 -1) dn and ( 34)

W = klmn c klmn I a k 1 -3 a k I b l 2 -3 b l I cm 4 -1 I dn 5 -1 , ( 35
)
where c klmn , a k , b l , c m and d n are material parameters. In the same way, other power law constitutive equations were proposed and are listed in Table 5. Additional forms can be found in the polyconvex strain energies listed in Table 1. These models represent different forms that link different invariants.

Exponential Development

A key property of the constitutive equation for soft tissues is the inclusion of an important strain hardening. This is easily obtained by means of an exponential function of the I 4 invariant. This approach is largely used in the literature, the first models were proposed in the 1990s. In the beginning, two fibre directions were introduced to represent the mechanical behaviour of arteries [104]. This was extended to four directions [START_REF] Baek | Theory of small on large: potential utility in computations of fluid-solid interactions in arteries[END_REF]159] and to n directions [START_REF] Gasser | A rate-independent elastoplastic constitutive model for biological fiberreinforcedcomposites at finite strains: continuum basis, algorithmic formulationand finite element implementation[END_REF] and used for example with 8 directions for cerebral aneurysms [276]. These models may be used to model the behaviour of a complex tissue such as in different areas of a soft tissue (as for example the different layers of an artery) [START_REF] Balzani | On the mechanical modeling of anisotropic biological soft tissue and iterative parallel solution strategies[END_REF]. Various formulations are listed in Table 6.

In order to take into account the ratio of isotropic to anisotropic parts of a heterogeneous material, a weighting factor has been introduced based on the contributions of I 1 and I 4 [107]. This represents a measure of dispersion in the fibre orientation. This model leads to

W aniso = c 1 2c 2 exp c 2 (I 4 -1) 2 -1 Weiss et al. [268] W aniso = c 1 exp(I 4 -1) 2 -(I 4 -1) 2 -1 Peña et al. [196] W aniso = c 1 c 2 exp c 2 (I 4 -1) -c 2 (I 4 -1) -1 I 1 , I 4 forms Holzapfel et al. [105] W aniso = c 1 (I 4 -1) exp c 2 (I 4 -1) 2 Gasser et al. [83] W = c 1 2c 2 [exp c 2 κI 1 + (1 -3κ)I 4 -1 2 -1] Holzapfel et al. [107] W = c 1 2c 2 exp c 2 (1 -κ)(I 1 -3) 2 + κ(I 4 -1) 2 -1 May-Newman and Yin [161, 162] W = c 0 exp c 1 (I 1 -3) 2 + c 2 ( √ I 4 -1) 4 -1 Rubin and Bodner [221] W = c 1 2c 2 exp c 2 c 5 (I 1 -3) + c 3 c 4 ( √ I 4 -1) 2c 4 -1 Lin and Yin [148] W = c 1 exp c 2 (I 1 -3) 2 + c 3 (I 1 -3)(I 4 -1) + c 4 (I 4 -1) 2 -
W aniso = C 1 2C 2 exp C 2 (I 4 -1) 2 -1 + C 3 2C 4 exp C 4 (I 5 -1) 2 -1
the creation of different constitutive equations which are also listed in Table 6. Recently, a general form of an energy function was devised [197] in order to summarise a large number of constitutive equations:

W = γ aη exp η(I 1 -3) a -f 1 (I 1 , a) + c i bd i exp d i I (i) 4 -I 0 4 b -g I (i) 4 , I 0 4 , b . ( 36
)
The choice of the functions f 1 and g allows for the wide generalization of many different models. Also, γ , η, a, b, c i , d i and I 0 4 are material parameters, and I 0 4 represents the threshold to reach for the fibre to become active. Some authors [START_REF] Einstein | Inverse parameter fitting of biological tissues: a response surface approach[END_REF][START_REF] Freed | Invariant formulation for dispersed transverse isotropy in aortic heart valves[END_REF] have proposed in a way similar as to what is done in the case of isotropy [START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF] a constitutive equation for the stress, the energy being obtained by integration: Horgan and Saccomandi [111] W ,4 = -

σ (λ) = A exp B λ 2 -1 2 -1 , ( 37
)
W = 2 3 I 4 c 2 1 + 2 c 1 √ I 4 - 3
Horgan and Saccomandi [111] W = -c 1 c 2 log 1 -

(I 5 -1) 2 c 2 Lurding et al. [153] W = c 1 (I 1 -3) + c 2 (I 1 -3)(I 4 -1) + c 3 (I 2 4 -I 5 ) + c 4 ( √ I 4 -1) 2 + c 5 ln(I 4 ) Chui et al. [153] W = c 1 ln(1 -T ) + c 5 (I 1 -3) 2 + c 6 (I 4 -1) 2 + c 7 (I 1 -3)(I 4 -1) with T = c 2 (I 1 -3) 2 + c 3 (I 4 -1) 2 + c 4 (I 1 -3)(I 4 -1)
Other invariants

Lu and Zhang

[151] W = c 2 exp c 1 ( √ I 4 -1) 2 + 1 2 c 3 (β 1 -1) + 1 2 c 4 (β 2 -2)
where A and B are material parameters. Even if this approach was initially developed and used for arteries [START_REF] Chagnon | An osmotically inflatable seal to treat endoleaks of type 1[END_REF]205,260], it is also often used for different living tissues, as for example human cornea [190], erythrocytes [130], the mitral valve [203], trachea [155,254], cornea [139,179], collagen [125], abdominal muscle [101].

Some Other Models

Other ideas have been developed for rubber like materials, as for example the Gent [START_REF] Gent | A new constitutive relation for rubber[END_REF] model which presents a large strain hardening with only two parameters. Its specific form gives it a particular interest for some tissues. This model was generalised to anisotropy in two ways [111]. Other different forms can be proposed with a logarithmic or a tangent function. A list of constitutive equations is given in Table 7. There are two ideas in these models. One is to describe the behaviour at moderate deformation. Thus, functions that provide for a weak slope are used; these models are principally used before the activation of muscles, i.e., when the material is very soft. When the material becomes stiffer, a function that models a large strain hardening is necessary. In this way, different functions were introduced to capture very important changes of slopes.

Coupling Influence

Different coupling can be taken into account in the constitutive equation, for example, the shear between the fibres and the matrix, and the interaction between the fibres.

Fibre Shear. In this case, the soft tissue is considered as a composite material, the strain energy is decomposed into three additive parts W = W m + W f + W f m [199], where the three terms are the strain energy of the matrix, of the fibres and of the interactions between fibres and matrix, respectively. Moreover, the deformation gradient of the fibres F can be decomposed into a uniaxial extension tensor F f and a shear deformation F s , as F = F s F f [START_REF] Criscione | Physically based strain invariants set for materials exhibiting transversely isotropic behavior[END_REF]. The decomposition of the strain energy function into different parts allows, for different loading states, the consideration of constitutive equations which are specific for the strain endured by the fibre, the matrix and the interface. This leads to the construction of different function forms [START_REF] Guo | Large deformation response of a hyperelastic fibre reinforced composite: theoretical model and numerical validation[END_REF]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ W m = 1 2 c 1 f (I 4 )(I 1 -3), W f = c 1 g 1 (I 4 ) I 5 -I 2 4 I 4 , W f m = c 1 g 2 (I 4 ) I 1 - I 5 + 2 √ I 4 I 4 . ( 38
)
Another basic form also has been proposed [199]:

W f m = g 2 (I 4 ) I 4 I 3 (I 5 -I 1 I 4 + I 2 ) -1 2 , ( 39
)
where f , g 1 and g 2 are functions to define and c 1 is a material parameter. The first function corresponds a generalisation of the neo-Hookean model [START_REF] Guo | Mechanical response of neo-Hookean fiber reinforced incompressible nonlinearly elastic solids[END_REF]. Few functions for f , g 1 and g 2 have so far been proposed, the first being based on exponential functions [START_REF] Caner | Hyperelastic anisotropic microplane constitutive model for annulus fibrosus[END_REF][START_REF] Guo | A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus[END_REF][START_REF] Guo | Large deformation response of a hyperelastic fibre reinforced composite: theoretical model and numerical validation[END_REF].

Interaction Between Fibres. Few models are proposed to take into account the influence of the coupling between different fibre directions. Different techniques can be used.

In order to take into account different directions and to not limit the problem to one direction fibre, it is also possible to couple invariants from different directions [228], the following invariant expression has been proposed: (1) 4

α 2 I (1)2 4 + 2λ(1 -α)I
I (2) 4 + (1 -α) 2 I (2)2 4 with α ∈ [0, 1]. ( 40
)
α represents a material parameter. This expression takes into account the deformation in two directions with only one invariant. Nevertheless, this has not yet been used in constitutive equations.

Instead of employing an additive decomposition of the strain energy to account for the different directions, a function that represents a coupling between the invariants of different directions [102] can be used [242]:

W = c 1 c 2 exp c 2 I (1) 4 + I (2) 4 -2 -c 2 I (1)
4 + I (2) A generalised weighted expression of the constitutive equation also has been developed [START_REF] Ehret | A polyconvex hyperelastic model for fiber-reinforced materials in application to soft tissues[END_REF]121]:

W = 1 4 r μ r 1 α r exp α r i γ i I (i) 4 -1 -1 + 1 β r exp β r i γ i J (i) 5 -1 -1 . ( 42
)
Even if the model was developed for pneumatic membranes, such representations that have proposed multiplicative terms between the I 4 invariants of each direction instead of an additive decomposition can be used for soft tissues [215]: (1,2) c I (1) 4 -

W = c (1) 1 I (1) 4 -1 β 1 + c (1) 2 I (1) 5 -1 β 2 + c (2) 1 I (2) 4 -1 γ 1 + c (2) 2 I (2) 5 -1 γ 2 + c (1) c (I 1 -3) δ 1 I (1) 4 -1 δ 1 + c (2) c (I 1 -3) δ 2 I (2) 4 -1 δ 2 + c
1 η I (2) 4 -1 η . ( 43
)
This strain energy introduces coupling between the different directions, but the additive decomposition of the constitutive equation allows one to fit separately the different parameters c (j ) i , δ i and β i .

Use of I 8 and I 9 . As proposed in the first part of this paper, coupling terms including I (i,j ) 8

and I (i,j ) 9 can be used. Thus such terms have been added to the strain energy in order to model esophageal tissues [177]:

W = c 1 c 3 exp c 3 (I 1 -3) + c 2 c 5 exp c 5 (I 2 -3) + c 4 c 2 7 exp c 7 I (1) 4 -1 -c 7 I (1) 4 -1 -1 + c 6 c 2 8 exp c 8 I (2) 4 -1 -c 8 I (2) 4 -1 -1 + c 9 I (1,2) 8 -I (1,2) 9 2 , ( 44
)
where c i with i = 1...9 are material parameters. For annulus fibrous tissues, the influence of the interaction between the layers has been modelled [178] with an energy term taking into account I (1) 4 , I (2) 4 and I (1,2) 8 :

W = c 1 2c 2 exp c 2 I (1,2)
8 (I (1) 4 I (2) 4 I (1,2) 9

) 1/2 -I (1,2) 9 2 -1 . ( 45
)
A similar form of exponential model (cf. Table 6) has been proposed to include the effect of I 8 [START_REF] Göktepe | Computational modeling of passive myocardium[END_REF]:

W = c 1 c 2 exp c 2 (I (1,2) 8 ) 2 I (1,2) 9 -1 . ( 46
)
These models are not often employed, but there exist some for composite materials that can be used [200,214]. In comparison with other models, these approaches take into account the shear strain in the material whereas the first models couple the deformations of the different fibres.

Statistical Approaches

In this part, some statistical approaches that tend to encompass the physics of soft tissues physics are detailed. They come from the study of the collagen network and use a change of scale method [START_REF] Chen | Nonlinear micromechanics of soft tissues[END_REF]181]. A collagen molecule is defined by its length, its stiffness and its helical structure. Some studies are motivated by approaches developed for rubber like material [START_REF] Beatty | An average-stretch full-network model for rubber elasticity[END_REF]73,129]. Unlike polymer chains in rubber which are uncorrelated in nature, collagen chains in biological tissues are classified as correlated chains from a statistical point of view. Rubber chains resemble a random walk whereas biological chains often present privileged oriented directions. It this way, different theories are considered to represent the chains, as for example wormlike chains with a slight varying curvature [132], or sinusoidal, zig-zag or circular helix representations [START_REF] Freed | Invariant formulation for dispersed transverse isotropy in aortic heart valves[END_REF]126,140]. Nevertheless, to develop models which rest on statistical approaches, some hypotheses are needed. A distribution function f of the orientation of the fibres is used to represent the material. The unit vector a 0 oriented in the direction of a certain amount of fibres having a spatial orientation distribution f is defined in terms of two spherical angles, denoted as φ and ψ:

a 0 = sin φ cos ψe 1 + sin φ sin ψe 2 + sin φe 3 , ( 47
)
with φ ∈ [0, π] and ψ ∈ [0, 2π] and e i is the usual rectangular Cartesian basis. The distribution function is required to satisfy some elementary properties [191]. By symmetry requirements f (a 0 ) = f (-a 0 ). The quantity f (a 0 ) sin φdφdψ represents the number of fibres with an orientation in the range [(φ, φ + dφ), (ψ, ψ + dψ)]. By considering the unit sphere S around a material point, the following property is deduced:

1 4π S f (a 0 )dS = 1 4π π 0 2π 0 f (a 0 ) sin φdφdψ = 1. (48)
A constant distribution leads to isotropy [START_REF] Ateshian | Anisotropy of fibrous tissues in relation to the distribution of tensed and buckled fibers[END_REF].

The strain energy of the soft tissue can then be deduced by integration of the elementary fibre energy in each direction w(I 4 (a 0 )) by:

W = 1 4π S f (a 0 )w I 4 (a 0 ) dS. ( 49
)
Finally, the stress is determined by derivation:

S = 1 2π S f (a 0 ) ∂w(I 4 (a 0 )) ∂C dS. ( 50
)
The evaluation of the stress depends on different parameters: the distribution function and the energy of a single fibre. Different considerations have been proposed in the literature. For the distribution function, the principal propositions are: beta distribution [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Cacho | A constitutive model for fibrous tissue sconsidering collagen fiber crimp[END_REF]218,225], log-logistic distribution [277], Gaussian distribution [START_REF] Billar | Biaxial mechanical properties of the native and glutaraldehyde-treated aortic valve cusp: Part II-a structural constitutive model[END_REF][START_REF] Chen | The structure and mechanical properties of the mitral valve leaflet-strut chordae transition zone[END_REF][START_REF] Driessen | A structural constitutive model for collagenous cardiovascular tissues incorporating the angular fiber distribution[END_REF]141,223,272], von Mises distribution [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF][START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF][START_REF] Girard | Peripapillary and posterior scleral mechanics-Part I: development of an anisotropic hyperelastic constitutive model[END_REF]191,211,263] or the Bingham distribution [START_REF] Alastrué | On the use of the bingham statistical distribution in microsphere-based constitutive models for arterial tissue[END_REF]. The forms of the distribution are listed in Table 8. The choice of the functions is also a key point. Different functions can be chosen to describe the mechanical behaviour of a collagen fibre; the simple linear behaviour [START_REF] Ateshian | Anisotropy of fibrous tissues in relation to the distribution of tensed and buckled fibers[END_REF], or the phenomenological laws of the exponential Fung type [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF][START_REF] Billar | Biaxial mechanical properties of the native and glutaraldehyde-treated aortic valve cusp: Part II-a structural constitutive model[END_REF]125,211,225,263] or a logarithmic function [277] or a polynomial function [START_REF] Flynn | An anisotropic discrete fibre model based on a generalised strain invariant with application to soft biological tissues[END_REF]248], other functions [119,207] or worm-like chain forms [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF][START_REF] Alastrué | On the use of the bingham statistical distribution in microsphere-based constitutive models for arterial tissue[END_REF][START_REF] Bischoff | A microstructurally based orthotropic hyperelastic constitutive law[END_REF][START_REF] Bischoff | Orthotropic hyperelasticity in terms of an arbitrary molecular chain model[END_REF][START_REF] Bustamante | Ten years of tension: single-molecule DNA mechanics[END_REF][START_REF] Garikipati | A continuum treatment of growth in biological tissue: the coupling of mass transport and mechanics[END_REF]135,136, 216,218] which are a particularisation of the eight-chain model [START_REF] Arruda | A three dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF] to the transversely isotropic case. For some models, a parameter should be introduced in the fibre concentration factor to control collagen fibre alignment along a preferred orientation [START_REF] Girard | Peripapillary and posterior scleral mechanics-Part I: development of an anisotropic hyperelastic constitutive model[END_REF]. The different constitutive equations are listed in Table 9. The reader can refer to [START_REF] Bischoff | Continuous versus discrete (invariant) representations of fibruous structure for modeling non-linear anisotropic soft tissue behavior[END_REF] to determine which strain energy is used for each tissue.

distribution β(η, γ ) = Γ (η)Γ (γ ) Γ (η+γ ) with Γ (x) = ∞ 0 t x-1 exp(-t)dt Log-logistic distribution f (ε) = k b (ε-ε 0 /b) k-1 [1+(ε-ε 0 /b) k ] 2 with ε = √ I 4 -1 Gaussian distribution f (φ) = 1 σ √ 2π exp --(φ-M) 2 2σ 2 Normalized von Mises distribution f (φ) = 4 I b 2π exp b(cos 2φ + 1) with I = 2 √ π √ 2b 0 exp -t 2 dt Bingham distribution f (r, A) = K(A) -1 exp r T • Ar A is a symmetric matrix, r a vector and K(A) a normalized constant
w = k 1 k 2 exp k 2 (I 4 -1) 2 -1 Logarithmic function w = c ε -log(ε + 1) for ε > 0 with ε = √ I 4 -1 Polynomial function w = 1 2 K γ + M m=2 γ m m γ γ m m Worm-like chain w = nkθ L 4A 2 r 2 i L 2 + 1 1-r i /L - r i L - ln(I 2 4 r 2 0 ) 4r 0 L 4 r 0 L + 1 [1-r 0 /L] 2 -1 -W r with r i = √ I 4 r 0 and W r = 2 r 2 0 L 2 + 1 1-r 0 /L - r 0 L 164,
The main difficulty of the different constitutive equations is that they need a numerical integration that is always time consuming [START_REF] Bischoff | Continuous versus discrete (invariant) representations of fibruous structure for modeling non-linear anisotropic soft tissue behavior[END_REF]211]. The integration of the fibre contribution is mainly realised over a referential unit sphere [134,172]. Some prefer to use a finite number of directions, the constitutive equation is thus modified as follows:

1 4π S (•)dS = m i=1 w i (•) i . ( 51
)
Different choices exist, as the 42 directions of Bazant and Oh [START_REF] Bažant | Efficient numerical integration on the surface of a sphere[END_REF] and by Menzel [164], or the 184 directions of Alastrue et al. [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF], for example. The only different approach to those mentioned above is that proposed by [START_REF] Flynn | An anisotropic discrete fibre model based on a generalised strain invariant with application to soft biological tissues[END_REF] who used only six initial directions without employing an integration. Even if statistical approaches have more complex equations than phenomenological ones, some of these models have been implemented in finite element codes [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF][START_REF] Garikipati | A continuum treatment of growth in biological tissue: the coupling of mass transport and mechanics[END_REF]135,136,164,275].

Discussion

The main difficulty is not to find the better constitutive equation but to have suitable experimental data. In fact, the difficulty is often that there is a large dispersion in experimental data due to the dispersion between the different specimens. Moreover, it is often difficult to impose different loading conditions on similar specimens. Thus the errors are often large and the number of loading conditions is often limited. As a consequence, one can wonder if the key point is to obtain the best fit for a very specific experimental database, or if the most important point is to represent globally the mechanical behaviour keeping in mind the physics of soft tissues.

As it was shown in the previous paragraphs, the number of constitutive equations that can be used to describe soft tissues non-linear elasticity is very impressive. Moreover, there exist other approaches, not presented in this paper, which involve a new class of elastic solids with implicit elasticity [212] that can also describe the strain limiting characteristics of soft tissues [START_REF] Freed | An implicit elastic theory for lung parenchyma[END_REF]. These theories are also elastic as they do not dissipate energies even though they are written in terms of strain rate and stress rate. But in this paper, we only focus on hyperelastic energy functions. These functions are expressed in terms of strain tensor components or strain invariants. The main difference between the two approaches discussed here is that the invariants formulation permits one to split the energy function into additive isotropic and anisotropic parts, even if, some constitutive equations written in invariants also link these two parts. The first constitutive equations introduced for soft tissues were isotropic. Although, for some applications, an isotropic constitutive equation is used to describe the mechanical behaviour for different soft tissues, the use of such simplified models is, in many cases, misleading and inappropriate as most soft tissues have a fibre structure that must be taken into account. To represent this structure, many constitutive equations are based on privileged directions that correspond to physical fibre orientations. In the modelling, characteristic directions are defined and they are represented by an angle that defines the orientation of the fibre compared to a specific direction. This angle can be considered as a parameter that is used to fit as well as possible the experimental data. Thus, the model is not used to mimic the physical soft tissue but it is used as a phenomenological equation to describe properly experimental data. This is not, in our opinion, a good choice, and it may mean that the energy function is not well chosen. The angle between the fibres should not be an adjustable parameter but must be imposed by the soft tissue structure.

An important issue in modelling concerns the stretching resistance of fibres. Many authors consider that the fibre must reach a threshold before opposing a stress. In this way, a threshold parameter can be introduced in all the suitable constitutive equations presented in this review. For the phenomenological model, it consist in replacing (I 4 -1) by (I 4 -I 0 4 ), or ( √ I 4 -1) by ( √ I 4 -I 0 4 ) in the constitutive equations. I 0 4 corresponds to the needed deformation to generate stress, see for example [START_REF] Calvo | On modelling damage process in vaginal tissue[END_REF]107,197,219]. The advantage of such approaches is that there is a material parameter that controls the beginning of material stiffening. Nevertheless, a main difficulty is that it strongly depends on the zero state of the experimental data. Moreover, this zero state is often different between post-mortem and in-vivo specimens, and can depend on the experimenter.

Anisotropic strain energy functions are difficult to fit, as it is difficult to separate the contribution between the matrix and the fibres, and to distinguish the different parts of the strain energy. Nevertheless, some strategies based on dissociating isotropic and anisotropic parts can be used [START_REF] Harb | A new parameter identification method of soft biological tissue combining genetic algorithm with analytical optimization[END_REF]. To avoid such representations, physical approaches attempt to represent the repartition of fibres in space, but two difficulties must be considered; the knowledge of the distribution function of the fibres in space and the mechanical properties of a single fibre. The choice of the best strain energy function is always a difficult point in the modelling process. A summary of the constitutive equations is presented in Fig. 1. In practise, the invariants I 2 and I 5 are often neglected. Their contribution is always difficult to determine [115] but it can be useful [START_REF] Feng | Measurements of mechanical anisotropy in brain tissue and implications for transversely isotropic material models of white matter[END_REF]. Moreover, these invariants are not independent from I 1 and I 4 in uniaxial loading tests. In this case, it is important to have also biaxial loadings to fit constitutive equations [224]. Moreover, in vivo experimental data [233] would be a benefit to obtain a good experimental fit, but there is little such data in the literature as compared to post-mortem experimental data. The constitutive equation choice will depend on the particular soft tissues under study and the conclusions will strongly depend on the experimental data that is chosen. Nevertheless, some comparisons between anisotropic strain energies have been realised in particular cases, see, for example, [START_REF] Carboni | Passive mechanical properties of porcine left circumflex artery and its mathematical description[END_REF][START_REF] Galle | A transversely isotropic constitutive model of excised guinea pig spinal cord white matter[END_REF]104,116,117,265].

In practice, a strategic point is the choice of a constitutive equation that is implemented in a finite element code to describe loading conditions that are very far from uniaxial or biaxial loadings. In this case, it is important to choose a constitutive equation that can be fitted with few experimental data that do not simulate non-physical response for any loading. Generally, it is better to limit the number of invariants and material parameters. Moreover, the simplest functions are often the best as they stand the least probability of creating nonphysical responses even if their fitting is not the best.

Conclusion

This paper has listed many different constitutive equations that have been developed for soft tissues. The number of constitutive equations to represent the contribution due to hyperelasticity is extensive due to the number of soft tissues and the experimental data dispersion. The paper has listed first, isotropic constitutive equations, and next anisotropic ones, and these were classed in different categories; those written with strain tensor components, those written in terms of the invariants, and those based on statistical modelling. Despite all the difficulties encountered in the modelling of the isotropic or anisotropic hyperelastic behaviour of soft tissue, these constitutive equations must be considered as only the basis of a more complex constitutive equation. Generalized equations should take into account other phenomena such as the activation of muscle [START_REF] Calvo | Passive non linear elastic behaviour of skeletal muscle: experimental results and model formulation[END_REF]158,188,189,250] or the viscoelasticity of the tissues [START_REF] Bischoff | A rheological network model for the continuum anisotropic and viscoelastic behavior of soft tissue[END_REF][START_REF] Haslach | Nonlinear viscoelastic, thermodynamically consistent, models for biological soft tissue[END_REF]105,147,208] or stress softening [154,195], for example. Nevertheless, the hyperelasticity representation should remain as the starting point in a modelling program and should be described as well as possible before introducing other effects.

2 exp c 2 (I 4 - 1 )

 2241 c 2 (I 4 -1) -1Holzapfel et al. [104]

1 2 + c 5 (I 1 - 3 ) 1 + c 5 (

 1251315 Doyle et al. [64] W = c 1 exp c 2 (I 1 -3) 2 + c 3 (I 1 -3)(I 4 -1) + c 4 (I 4 -1) + c 6 (I 4 -1) -1 Fung et al. [78] W = c 1 exp c 2 (I 1 -3) 2 + c 3 (I 1 -3)(I 4 -1) + c 4 (I 4 -1) 2 -I 1 -3) 2 + c 6 (I 4 -1) 2 + c 7 (I 1 -3)(I 4 -1) I 4 , I 5 forms Masson et al. [160]

c 1 1 ) 4 +

 114 +c 2 (I 4 -1) c 3 +c 4 (I 4 -Ruter and Stein [222] W = c 2 cosh(I 4 -1) -1 Horgan and Saccomandi [111] W = -c 2 c 3 I 4 -1 + c 3 ln(1c 6 ln(I 4 ) Calvo et al. [38] W = 2c 5 √ I 4 + c 6 ln(I 4 ) + c 7Demirkoparan et al.[START_REF] Demirkoparan | Swelling of an internally pressurized nonlinearly elastic tube with fiber reinforcing[END_REF][START_REF] Demirkoparan | On dissolution and reassembly of filamentary reinforcing networks in hyperelastic materials[END_REF]

Fig. 1

 1 Fig. 1 Organisation of the constitutive equations in the paper

Table 2

 2 Principal isotropic hyperelastic constitutive equations developed for soft tissues, where c 1 , c 2 , c 3 and c 4 are material parameters. (*) The model is known as the generalised neo-Hookean model. (**) As pointed out by[109] is frequently mistakenly referred toDelfino et al. [56]

	Polynomial models
	Raghavan and Vorp [210]

Table 3

 3 Anisotropic constitutive equations written with strain tensors components, where A ij kl with i, j, k, l = 1...3, b i , with i = 1...12, a ij , b ij , c ij with i, j = 1...3 and c are material parameters

	Generic Fung functions W = C 2 (exp Q -1)
	Tong and Fung [252]

  E 11 + b 2 E 22 + b 3 E 33 + b 4 E 11 E 22 + b 5 E 11 E 33 + b 6 E 22 E 33 E 12 E 21 + b 10 E 3 11 + b 11 E 3 22 + b 12 E 2 11 E 22 + b 13 E 11 E 2

			2 zz
		+ b 7 E 3 zz
	Rajagopal et al. [213] Tong and Fung [252]	W = b 1 + b 7 E 2 11 + b 8 E 2 22 + b 9 E 2 33 + b 10 E 2 12 + b 11 E 2 13 + b 12 E 2 23 W = 1 2 b 1 E 2 11 + b 2 E 2 22 + b 3 E 2 12 + E 2 21 + 2b 4 E 12 E 21
		+	b 5 2 exp b 6 E 2 11 + b 7 E 2 22 + b 8 E 2 12 + E 2 21 + 2b 9 22 -1

  For example, the function defined in [175] requires a limit for each component. As a consequence, the domain limits of the function are well established. The question is different for[249] as the function is written in terms of the sum of the components in a logarithmic form and the function can be undefined[104].

	Nash and Hunter [175]	W = c 11	E 2 11 |a 11 -E 11 | b 11	+ c 22	E 2 22 |a 22 -E 22 | b 22	+ c 33	E 2 33 |a 33 -E 33 | b 33
		+ c 12	E 2 12 |a 12 -E 12 | b 12	+ c 13	E 2 13 |a 13 -E 13 | b 13	+ c 23	E 2

Table 4

 4 Some constitutive equations based on truncations of the series developments where c i with i = 1...6 are material parameters

	I 4 forms
	Triantafyllidis and Abeyaratne [257]	W aniso = c 2 (I 4 -1) 2
	Peng et al. [199]

Table 5

 5 Model based on a power development, where k i , i = 1..6 are material parameters

	Power developments
	Ghaemi et al. [85]	W aniso = C I	k 1 /2 4	-1	k 2

Table 6

 6 List of exponential constitutive equations, where c 1 , c 2 , c 3 , c 4 , c 5 and κ are material parameters

Table 7

 7 Other models written in invariants, where c i with i = 1...7 are material parameters 1 +c 2 (I 4 -1) c 3 +c 4 (I 4 -1)+c 5 (I 4 -1) 2 +c 6 (I 5 -2I 4 +1)

	General forms
	Horgan and Saccomandi [111]	W ,4 = -

c

Table 8

 8 Some distribution functions used in statistical approaches, where ε 0 , b, σ , M and I are statistical parameters

	Distribution functions
	Beta

Table 9

 9 Some fibre functions used in statistical approaches where k 1 , k 2 , r 0 , L, K, W r , γ m and m are material parameters

	Energy functions
	Holzapfel et al. function

|a 23 -E 23 | b 23Moreover, the parameters of these materials are often difficult to fit as they have no physical meaning. For example, the strain energy of [259] is discussed in [104] and is not convex, this can also be the case for Fung functions if the parameters are not well chosen[104]. The limitations in material parameters are discussed in[START_REF] Federico | An energetic approach to the analysis of anisotropic hyperelastic materials[END_REF] 269] with respect to polyconvexity. In this way, developments have been made to ensure polyconvexity with a physical meaning of the material response[247]. Other conditions also must be respected for viable functions.

+ 2c 2 -1 . (41)

Acknowledgements

The authors thank Prof. Roger Fosdick for his valuable comments. This work is supported by the French National Research Agency Program ANR-12-BS09-0008-01 SAMBA (Silicone Architectured Membranes for Biomedical Applications).