{'fraction_non_alphanumeric': 0.03186658826603946, 'fraction_numerical': 0.012838102704821639, 'mean_word_length': 4.855489671002295, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 1, 'zh_pornsignals.json': 0}}
An in silico-based review on anisotropic hyperelastic constitutive models for soft biological tissues

We review nine invariant and dispersion-type anisotropic hyperelastic constitutive models for soft biological tissues based on their fitting performance to experimental data from three different human tissues. For this, we used a hybrid multi-objective optimization procedure. A genetic algorithm is devised to generate the initial guesses followed by a gradient-based search algorithm. The constitutive models are then fit to a set of uniaxial and biaxial tension experiments conducted on tissues with differing fiber orientations. For the in silico investigation, experiments conducted on aneurysmatic abdominal aorta, linea alba, and rectus sheath tissues are utilized. Accordingly, the models are ranked with respect to an objec tive normalized quality of fit metric. Finally, a detailed discussion is carried out on the fitting performance of each model. This work provides a valuable quantitative comparison of various anisotropic hyperelastic models, the findings of which can aid those modeling the behavior of soft tissues in selecting the best constitutive model for their particular application.


Introduction
The physiological functions of many soft biological tissues are related to their material composition and mechanical properties. Examples of these functions include blood flow through the heart and blood vessels [12,59,68], peristalsis of the gastrointestinal (GI) tract [4,52], and even the transmission of sound within the ear cavity [90]. The mechanical properties of tissues can also be a predictor of their pathophysiology, such as the stiffness of arteries being a successful indicator of cardiovascular disease, or the elastic and viscoelastic parameters of skin being linked to the occurrence of autoimmune disorders [49,51]. To successfully model the physiological functions of soft tissues and to investigate any pathophysiology, the characteristics of the material must first be established. Constitutive modeling can be used to analytically capture the mechanical behavior of these biological tissues. Such an ansatz forms the basis for developing computational models which consider the geometry of the tissue or organ to study how shape contributes to its function. Baseline hyperelasticity is a precursor in biomechanics to the modeling of more advanced phenomena such as viscoelasticity [11], electroelasticity [19] and electro-viscoelasticity [9,10], to mention a few. In addition to this, computational modeling can be used to study fluid-structure interaction scenarios, such as urine traveling down the urinary tract [91], or air filling up the alveoli of the lungs [85]. However, before this is possible, the behavior of the solid structures must be determined and successfully simulated. A variety of applications exist for the computational models mentioned, including within the medical device industry, where the mathematical models can be used to aid in the design of devices that interact with the human body, such as stents and endoscopy devices [25,81]. The use of the models in a finite element implementation can give greater insight into the devices' mechanisms, highlight any potential risks for the patient, and conserve a range of resources such as time and manufacturing materials. Further to this, computational models of soft tissues can be used for surgical simulations; the training technique which allows surgeons to practice a procedure an infinite number of times without risking the health or lives of patients [1,32,61].
It has previously been well established that soft tissues respond elastically to a given force and past a certain point, i.e. under large strains, behave non-linearly [5,26,88]. For this, and when the time-dependent properties of the tissue are not considered, hyperelastic modeling, similar to that for polymers, is required [46,79]. For polymers, isotropic, hyperelastic models are commonly used as these substances often have a homogeneous material composition, unless specifically designed otherwise. Biological tissues, however, are generally observed through experimentation to be anisotropic. For fibrous tissues, this is often associated with different densities and/or orientations of the tissue's main constituent fibers, resulting in a differing response of the tissue depending on its loading direction. The anisotropy of biological tissues, hard or soft, can normally be linked with the physiological function of the organ when considering the mechanical characteristics of each different fiber. There are many types of tests to characterize the anisotropic hyperelastic behavior of biological tissues, and amongst the most common include uniaxial tension tests in two or more directions [62], biaxial tension tests which simultaneously load in two directions [29,47], pure shear tension tests in two directions [56], inflation-extension tests for tubular tissues [76] and compression tests in different directions [7]. Experiments conducted on soft biological tissues can be difficult due to their delicate nature in addition to practical limitations including the need for highly sensitive data acquisition equipment to capture their comparatively low internal forces. This makes the experimental studies conducted on soft tissues extremely valuable in the field of biomechanics. Such studies conducted using animal soft tissues are more prevalent than with human tissue due to the former being more readily available and having fewer associated ethical constraints. To date, some non-exhaustive examples of the experimental data available on animal soft tissues includes the porcine aorta [67], rodent and porcine brain [53], rodent vagina [47], porcine liver [57], ovine oesophagi [77], swine tongue [86], equine urethra [63], porcine skeletal muscle [80], rabbit skin [75] and porcine stomach [48].
Animal tissue may provide a good representation of how human tissue responds to certain loads or strains, but the data from such experiments cannot be used to accurately model human tissues, especially for use within medical applications. Despite being less prevalent than animal studies, there are also a wide variety of experimental studies on human soft tissues, with investigations being conducted both in vivo and ex vivo. In this review, we will focus purely on human data from ex vivo experiments as these allow for large deformations that clearly display the non-linear elastic behavior of the tissues. Experimental data for healthy human soft tissues can be found for the brain [8], spleen [50], skin [33], aorta [24,84], bile duct [31], oesophagus [21,22], GI tract [23], liver [83], cervical tissue [62], while also for diseased tissue such as aneurysmatic aortas [78] and diabetic foot plantar tissues [66]. As mentioned previously, and as can be witnessed from the timeindependent stress-strain relationships of these experimental studies, at large strains soft biological tissues exhibit a non-linear hyperelastic material response, therefore their behavior cannot be simply captured by two constants of linear elasticity. A strain energy density function (SEDF) must be used to model the hyperelastic response of a soft tissue to certain loads or deformations, and SEDFs exist in many different forms. There are a variety of approaches including the very popular classical models by Fung et al. [27] and Holzapfel et al. [36], as well as more recently developed models [14,15], and those which describe certain time-independent phenomena such as the Mullins effect [69]. Micromechanical models provide insight into the mechanical role of the different components that comprise a material. In these models, every parameter is linked to a physical phenomenon, with the aim to predict the behavior of the tissue based purely on its composition. Micromechanical dispersion models give further insight into the orientation and distribution of the tissues' main constituent fibers. The parameter in these models associated with the dispersion of the fibers should ideally be identified using microstructural data obtained from the tissue in question [38]. Phenomenological models, on the other hand, utilize material parameters that do not have a physical meaning, these, however, are often found to capture the experimental behavior very well and may be more relevant for a specific application or loading domain [27].
A number of comprehensive reviews have been carried out on the topic of hyperelastic constitutive models. Steinmann et al. [79], Hossain et al. [44,45] and Dal et al. [18] reviewed the range of hyperelastic models in literature. A review of the models used to simulate the behavior of anisotropic hyperelastic soft biological tissues was conducted by Chagnon et al. [13]. This review presents an extensive and detailed list of many prominent models used in the field. Wex et al. [89] provided a review of the most commonly used isotropic constitutive models for biological tissues. Mihai et al. [60] reviewed models used for the brain and fat tissues, and successfully compared their performance ability under different loading conditions. Furthermore, Rosen and Jiang [70] investigated current isotropic hyperelastic laws able to model tissue-mimicking materials characterized through shear wave elastography. Martins et al. [54] carried out a concise review on isotropic models subjected to uniaxial loading. Wang et al. [87] added a level of complexity and reviewed the state of modeling the compressibility of soft tissues. Holzapfel et al. [40] looked at purely dispersion-type models. Bhattarai et al. [6] reviewed the modeling of only one tissue, the large intestine. Recently, Cheng and Zhang [16] provided a framework for the derivation of hyperelastic models, both isotropic and anisotropic, in the context of biomaterials.
The aforementioned reviews exclude a one-to-one quantitative comparison of models used for soft tissues. In this work, we provide a derivation of invariant-and dispersion-type models, an outline of the most prominent constitutive models in literature, and an in silico comparison of the models' ability to simulate the anisotropic hyperelastic behavior of a variety of human tissues. The data for this is obtained from experimental studies of aortic aneurysms and different regions of the healthy abdominal wall. This work provides a detailed framework of how to utilize the constitutive models, while also providing an in-depth critical analysis of each models' performance ability according to different data sets and loading domains. The findings of this review can be used in research or industry to support and inform the selection of the best constitutive model for a specific soft tissue application.
We pursue our investigation based on the datasets proposed by Niestrawska et al. [64], Cooney et al. [17], and Martins et al. [55] for abdominal aortic aneurysm (AAA) tissue, the linea alba and the rectus sheath, respectively. AAA is a pathological condition of the abdominal aorta which results in local bulging, while the linea alba and rectus sheath are sections of the abdominal wall. One reason for selecting these tissues is the fact that they are often subjected to surgical treatment, therefore an accurate representation of their mechanical response leads to representative numerical analysis for surgical procedures and the design of proper medical devices. In addition, each tissue presents a unique stress-stretch response, thus the choice of these tissues allows the models to be assessed across different mechanical behaviors. The overall output of the paper is twofold. Firstly, the applicability of each constitutive model to a particular type of tissue will be assessed, specifically the aptness of various anisotropic fiber distribution functions to the three tissues mentioned above. Secondly, the fitting performance of various dispersion-type anisotropic constitutive models will be compared. In this study, we will • outline nine anisotropic hyperelastic models and derive the closed form stress expressions for the uniaxial tension (UT) and equibiaxial tension (ET) deformation states based on the classical incompressibility assumption, • implement the aforementioned constitutive models into a hybrid genetic-gradient search algorithm outlined in Dal et al. [18], • identify the material parameters of each model (i) with respect to an ET dataset for AAA tissue [64], (ii) a UT dataset for the linea alba [17], and (iii) a UT dataset for the rectus sheath [55], • sort all the models with respect to an objective quality of fit metric according to their fitting performance to the AAA tissue dataset [64], the linea alba dataset [17], and the rectus sheath dataset [55]. • The identified parameters and the quality of fit values of each constitutive model will be presented in tabular form. Also, the stress-strain curves for each constitutive model will be provided in separate graphs.
The paper is organized as follows: the mathematical preliminaries for the kinematics of incompressible anisotropic hyperelastic solids are presented in Section 2. In Section 3, nine anisotropic hyperelastic constitutive models for soft biological tissues are summarized. Section 4 outlines the experimental procedures for the characterization of the quasi-static hyperelastic response of the soft human tissues selected in this study. In Section 5, the parameter identification procedure and quality of fit metric are provided. Finally, a detailed discussion of the results is presented in Section 6.

Fundamentals of hyperelastic deformable solids
In this section, the kinematics, state variables and stress expressions for an anisotropic hyperelastic continuum will be introduced. The mathematical framework of invariant and dispersion-type anisotropic hyperelastic formulations will be briefly discussed.

Kinematics
Let the deformation map ϕ(X, t) represent the motion of a deformable solid body. It maps the reference/Lagrangian configuration X ∈ B 0 of material point onto the current/Eulerian configuration of material points x = ϕ t (X) at time t ∈ T ⊂ R + . The deformation gradient maps the unit tangent of the reference configuration onto the spatial counterpart in the current configuration where T X B 0 and T x B denote the tangent spaces in the reference and current configurations, respectively. Additionally, the co-tangent spaces in the reference and current manifolds are represented as T * X B 0 and T * x B, respectively. The normal map between the unit normals in the undeformed and deformed configurations is defined as where n is the unit normal vector in the reference configuration and N is the unit normal vector in the current configuration. In this sequence, let line, area, and volume elements in the Lagrangian configurations to be denoted as dX, dA, dV , respectively. The Eulerian counterparts of these elements are obtained through the deformation gradient F , its cofactor cof[F ] = det[F ]F −T and its Jacobian J := det[F ] mappings between the co-and contravariant objects in the Lagrangian and Eulerian manifolds, the undeformed B 0 and deformed B configurations are locally equipped with covariant G and current g metric tensors in the neighborhoods N X of X and N x of x, respectively. Here, the right Cauchy-Green tensor and the inverse of the left Cauchy-Green tensors are defined as the pull back of the current metric g and push-forward of the Lagrangian metric G, respectively. The inverse of the left Cauchy-Green tensor is commonly known as the Finger tensor and represented as b = c −1 .

Principal stretches and invariants
The principle of material objectivity and the principle of frame indifference require that the energy stored in the hyperelastic material be a function of either principal stretches or invariants of the right Cauchy-Green tensor. The spectral decomposition of right Cauchy-Green tensor and its cofactor reads where are the principal areal stretches, see Figure 2. Moreover, the three isotropic invariants of the right

Free energy function and the stress expressions
Hyperelastic materials are governed by a potential function, the Helmholtz free energy function that describes the energy stored as a result of mechanical deformation. Polymeric materials and soft biological tissues exhibit a distinct response to bulk deformation and shear-type deformations. In this context, it is common practice to split the deformation into dilatational and volume-preserving parts The corresponding deformation measure reads Based on (12) and (13), the free energy function is expressed as where U (J) andΨ(F , A 1 , A 2 ) represent the volumetric and isochoric response of the material, respectively. Further decomposition of the isochoric free energy function was suggested by Holzapfel and Weizsacker [42]. They proposed decomposition of the free energy function into isotropic and anisotropic parts, The wavy collagen fibers do not store energy under contractile deformations. Hence, the isotropic ground matrix is activeΨ iso under small deformations. On the other hand, the collagen fibers dominate the overall behavior at high stretch levels governed byΨ ani [36]. Following the arguments of Gültekin et al. [34,35], we employ the specific form of the free energy function leading toΨ Herein, theF is replaced with F for the sake of convenience. 1 Recall that, in the incompressible limit, the formulations (14) and (16) lead to identical results. In this investigation, the volumetric part of the free energy function U (J) will be excluded and the pressure term will be obtained from boundary conditions by enforcing exactly the incompressible deformation state under uniaxial and biaxial deformations. A canonical relation between the stresses and the free energy function can be established in the sense where S and τ are the second Piola-Kirchhoff and the Kirchhoff stresses, respectively. Thereafter, the Lagrangian moduli expression result from the second derivative ofΨ with respect to C,

Invariant-based anisotropic hyperelastic formulations
The Lagrangian and Eulerian stresses for invariant-based formulations of anisotropic elasticity result from the application of the chain rule 1 In the susequent treatment, the biological tissues will be assumed to behave perfectly incompressible.
The closed form expressions for the Lagrangian and Eulerian stress tensors (26) require the following derivatives Substitution of (21) into (26) 1 gives the invariant-based representation for the second Piola-Kirchhoff stress tensor The derivatives with respect to the spatial metric g follow from the push-forward operation of (21) Substitution of (23) into (26) 2 leads to the invariant-based representation for the second Kirchhoff stress expression with ψ i = ∂Ψ ∂I i and p = −J∂ J Ψ. The stress expressions given in (22) and (24) can be expressed in terms of index notation for the principal directions The Lagrangian moduli expression for invariant-based formulation can be written as Furthermore, moduli expression can be additively decomposed into isotropic and anisotropic parts as follows The directional stiffness (DS) can be defined as DS := n ⊗ n : C ani : n ⊗ n with n = cosαe 1 + sinαe 2 , where α ∈ [0, 2π].

Dispersion-type anisotropic formulations
The tissue is considered a fiber-reinforced composite with the fibers distributed within an isotropic matrix in the dispersion-type anisotropic formulation; the models developed within this framework accurately describe the effect of the structural arrangement of the fibers on the mechanical response. Dispersion-type anisotropic approaches utilize density distribution functions to represent the distributed fiber architecture of tissues. Let unit fiber direction r on a unit-sphere be given in the undeformed configuration. The fiber density in direction r is expressed with ρ(r). The unit fiber orientation vector can be expressed in spherical coordinates as follows in terms of Eulerian angles θ ∈ [0, π] and φ ∈ [0, 2π], see Figure 3. The Eulerian counterpart of r is derived as t = F r. In this part, we will summarize the two kinematic approaches for the modeling of dispersion-type anisotropy in soft biological tissues.
PSfrag replacements (i) Generalized structure tensor (GST) formulations: The generalized structure tensor is defined as and | S |= 4π for a unit sphere. Let h = F HF T be defined as the Eulerian counterpart of the generalized structure tensor. The Lagrangian and Eulerian stresses can be obtained by applying the chain rule where is the one-dimensional average anisotropic fiber strain akin to the Green-Lagrangian strain. The closed form expressions for the Lagrangian and Eulerian stress tensors require the following derivatives (23) and (34) into (31) and (32), and by rearranging the terms, finally gives the generalized structure tensor-based representation of stresses for an anisotropic hyperelastic solid with ψ f = ∂ E Ψ ani . In the principal directions, (35) can be written in terms of the principal stretches, The anisotropic part of the Lagrangian moduli expression can be written as (ii) Angular integration (AI) formulations: The total free energy of fibers is calculated by where S represents a unit sphere with | S |= 4π. Starting with the free energy (38), the Eulerian stress tensor follows by the application of the chain rule (39) gives the Kirchhoff stress expression with ψ f = ∂ λF ψ f ib . The general form of the Kirchhoff stress based on the angular integration formulation is given in (40). In the principal directions, (40) can be written in terms of the principal stretches,

Stresses under homogeneous deformations
In this section, we briefly outline the deformation state of anisotropic biological tissues subjected to uniaxial and equibiaxial stresses. The deformation gradient tensor of anisotropic materials cannot be expressed in a straightforward manner as can be done for isotropic materials since the amount of contraction in the transverse direction and the shear deformations depend on the anisotropic fiber structure. Therefore, the deformation gradient tensor should be obtained iteratively to satisfy equilibrium conditions, e.g. in case of uniaxial tension σ 22 = σ 33 = σ 12 = σ 13 = σ 23 = 0. To this end, we assume that the material is deformed in the principal directions corresponding to the principal axes of fiber orientations for AI-or GST-based models. For invariant formulations, the deformation axes are assumed to coincide with the symmetry axes for the two fiber family formulations, while in single fiber family formulations one of the principal deformation axes coincides with the fiber direction. Such an ansatz avoids a priori shear strains/stresses, leading to diagonal deformation and stress tensors. In this context, the deformation gradient and the nominal stress expression under homogeneous uniaxial, equibiaxial or pure shear deformations can be expressed as follows Uniaxial tension: For an incompressible anisotropic hyperelastic solid, the deformation and stress states under uniaxial tension are see Figure 4(a),(b). The first two invariants under uniaxial deformation read For symmetrically orthotropic fibers, the equalities I 4 = I 6 and I 5 = I 7 hold leading to where ϕ is the angle between the fiber and the symmetry axis e 1 . The components of the nominal stresses under uniaxial loading in the symmetry axes for the invariant-, generalized structure tensorand angular integration-based formulations read Equibiaxial tension: For an incompressible hyperelastic anisotropic solid, the deformation and stress states under equibiaxial tension are see Figure 4(c). The first two invariants under equibiaxial deformation read The invariants associated with the structural tensors are I 4 = I 6 = λ 2 cos 2 ϕ + λ 2 sin 2 ϕ and I 5 = I 7 = λ 4 cos 2 ϕ + λ 4 sin 2 ϕ .
The components of the nominal stresses under equibiaxial loading in the symmetry axes for invariant-, generalized structure tensor-and angular integration-based formulations read Inv :

Hyperelastic material models
In this section, we will review a common isotropic model and nine anisotropic hyperelastic models under two main categories: (i) strain invariant-based models, and (ii) fiber dispersion-based models. The free energy functions and necessary derivatives for the stress expressions are outlined. In order to emphasize the degree of anisotropy and the anisotropy distribution obtained from paramater identification process from each tissue is presented in terms of polar plots for each constitutive model investigated. The polar plots consider the density distribution (DD) function for the plane of interest along with the directional stiffness (DS) computed from equation (28). For the invariant based formulations, the polar plots of the directional stiffness obtained merely from the anisotropic part of the initial elastic moduli tensor is used for visualization of the plane of interest.

The neo-Hookean model
The invariant-based models assume perfect alignment of fibers embedded into an isotropic ground matrix. Models have been presented for two families of fibers symmetrically disposed.
The neo-Hookean model is the most fundamental model for hyperelastic constitutive models. Many researchers represent the isotropic ground matrix of soft tissues with the neo-Hookean model [2,3,36,41,88]. Therefore, it deserves a separate description. Based on Wall's treatment of elasticity of a molecular network, Treloar [82] proposed the following free energy function with µ = nk B θ. Herein, n is the volume specific chain density, k B is the Boltzmann constant and θ is the absolute temperature, and µ = nk B θ is the shear modulus. The nonzero derivative of ψ with respect to the invariants is

Newman-Yin (NY) model
Newman and Yin [58] proposed a exponential free energy form analogous to the one proposed by Fung et al. [27] to describe hyperelastic behavior of mitral valve tissue. They assumed that the material is transversely isotropic and that the free energy function depends on the first and the fourth invariants. They observed that for a constant I 4 , both ψ 1 and ψ 4 increase nonlinearly. Therefore, Newman and Yin proposed, Fig. 5 The polar plots of the directional stiffness (DS) for the NY model [58]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT datset for the rectus sheath [55].
Q is the quadratic function of the invariants. In order to have a strain energy increase with increasing I 4 , k 0 should be positive. Also, if the tissue is not able to support compressive load, k 2 should be positive. The original model was proposed for tissues with a single family of fibers. In this work, we have extended the formulation considering two families of fibers as given in (53) 2 by incorporating the latter term depending on I 6 . The derivatives of Ψ ani with respect to the invariants are The polar plots of the directional stiffness for the NY model obtained with the identified parameters for each tissue are depicted in Figure 5.

Holzapfel-Gasser-Ogden (HGO) model
Holzapfel et al. [36] proposed a constitutive model for arteries. They additively decomposed the isochoric free energy function into isotropic and anisotropic parts, They utilized the neo-Hookean model (51) for the isotropic part since collagen fibers are thought to not contribute to the mechanical behavior of the tissue at low pressures. The free energy stored by two families of collagen fibers is described as where k 1 > 0 is a stress-like parameter and k 2 > 0 is a dimensionless parameter. The collagen fibers do not support compressive stresses due to their wavy nature. Therefore, the Macauley brackets are utilized in order to filter out the tensile stretches. The derivatives of ψ i with respect to the i th invariants are The polar plots of the directional stiffness for the HGO model obtained with the identified parameters for each tissue are depicted in Figure 6. Fig. 6 The polar plots of the directional stiffness (DS) for the HGO model [36]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT datset for the rectus sheath [55].

Holzapfel-Sommer-Gasser-Regitnig (HSGR) model
Holzapfel et al. [41] proposed a free energy function of the form Ψ(I 1 , I 4 ) = Ψ iso (I 1 ) + Ψ ani (I 1 , I 4 ). The model uses the neo-Hookean model (51) as the isotropic part of the free energy function. The anisotropic part of the free energy function has a mixed representation frag replacements The polar plots of the directional stiffness (DS) for the HSGR model [41]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT datset for the rectus sheath [55].
where k 1 > 0 is a stress-like parameter and k 2 > is a dimensionless parameter. The measure of anisotropy parameter p ∈ [0, 1] interpolates between the contributions of the first I 1 and fourth invariants I 4 . It mimics the degree of fiber dispersion phenomenologically. The switch function sgn I 4 − 1 enforces the tension-only condition and is activated for I 4 − 1 > 0. The model recovers the HGO model for p = 1. The fiber related terms drop and the model reduces to an exponential isotropic model for p = 0. The nonzero derivatives of the free energy function are  Fig. 8 The polar plots of the directional stiffness (DS) for the OS model [65]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT datset for the rectus sheath [55].

Ogden-Saccomandi (OS) model
Ogden and Saccomandi [65], following the work of Horgan and Saccomandi [43], proposed a logarithmic constitutive law for arterial tissue with two fiber families in which the fiber extension is limited. They additively decomposed the free energy function into isotropic and anisotropic parts, To model the isotropic behavior of the tissue, they adapted the well-known rubber elasticity model of Gent [30]. The free energy function of the Gent model reads where µ is the shear modulus and J m is the parameter that controls the chain extensibility limit for the matrix material. The deformation limit for the first invariant is I 1 < 3 + J m and the stresses tend to infinity asymptotically at this limit. As J m → ∞, the isotropic part of the free energy function (62) where k 1 is a stress-like parameter, and J f is the limiting parameter of extensibility of collagen fibers. The constraint I α < J f + 1 , α = 4, 6 .
(64) confines the stretches beyond the fiber extensibility limit. The derivatives of the free energy function are The polar plots of the directional stiffness for the OS model obtained with the identified parameters for each tissue are depicted in Figure 8.

Dispersion-type anisotropic constitutive models
It has been determined that not only do the mean orientation of fibers within soft tissues affect their mechanical properties, but also the amount by which they are dispersed around this mean [40]. For instance, if the majority of fibers are orientated along the direction of the mean, the behavior will be very different to that if the fibers are distributed in a cone-like span around the mean orientation [28]. For this reason, a number of models have been developed to incorporate this dispersion into the constitutive law describing the stress-strain relation of the tissue. These fiber dispersion-based models utilize a probability distribution function to model the histological structure of tissues. The models which have an angular integration approach and a generalized structure tensor approach are outlined in this section.

GST-based Gasser-Ogden-Holzapfel (GOH) model
Gasser et al. [28] assumed that the tissue's fibers are distributed rotationally symmetric around a mean fiber orientation direction M . They utilized a planar π-periodic von-Mises distribution as a fiber density distribution around M . The von-Mises distribution function is a one-dimensional probability distribution which is a function of Θ and concentration parameter b. The standard von-Mises distribution function is where I 0 (b) is a modified Bessel function of the first kind of order zero. Applying the normalization condition to (66) gives the relation The normalized von-Mises distribution is By inserting the von-Mises type density distribution function (68) into (30), the generalized structure tensor can be written as is the fiber dispersion parameter. A one-to-one relation exists between the dispersion parameter κ and the concentration parameter b, see Figure 9. Hence, κ ∈ [0, 1/3] enters the constitutive model as an additional material parameter responsible for the degree of dispersion. The lower limit κ = 0 recovers the invariant-based anisotropy and the upper limit κ = 1/3 leads to an isotropic constitutive response. The polar plots of the density distribution and directional stiffness for the PSfrag replacements Fig. 10 The polar plots of the density distribution (DD) and the directional stiffnesses (DS) for the GOH model [28]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT dataset for the rectus sheath [55]. GOH model obtained with the identified parameters for each tissue are depicted in Figure 10. Gasser et al. [28] can be considered as the GST counterpart of the GHO model [36]. The anisotropic part of the free energy function reads where E i = H i : C − 1 replaces the fourth and sixth invariants in the GHO model. The classical neo-Hookean model (51) is utilized for the isotropic part of the free energy function. The derivatives of the free energy function are

GST-based Holzapfel-Niestrawska-Ogden-Reinisch-Schriefl (HNORS) model
Holzapfel et al. [37] take into account both the in-and out-of-plane dispersion of fibers based on the observations of Schriefl et al. [71][72][73][74] where they recorded that the fibers are dispersed both in-plane and out-of-plane. Their observations reveal that no correlation exists between in-and outof-plane dispersions. Based on these arguments, the probability density function is multiplicatively exp(x cos α)dα (73) where a is the concentration parameter and I 0 (a) is the modified Bessel function of the first kind of order zero. The out-of-plane distribution is in the form The measures of dispersion in the in-plane and the out-of-plane directions read The structure tensor H has the form where M i is the in-plane mean fiber direction, whereas M n is the out-of-plane vector. The Eulerian counterpart of the generalized structure tensor reads where m i = F M i and m n = F M n , respectively. The lower and upper bounds for the in-plane and out-of-plane dispersion parameters are, respectively, 0 ≤ κ op ≤ 1/2 and 0 ≤ κ ip ≤ 1 .
The polar plots of the density distribution and directional stiffness for the HNORS model obtained with the identified parameters for each tissue are depicted in Figure 12. The free energy functions are identical to the GOH model [28] as depicted in (70).

AI-based Alastrué -Martinez-Doblaré -Menzel (AMDM) model
The angular integration-based anisotropic model of Alastrué et al. [2] takes into account rotationally symmetric fiber dispersion based on the micro-sphere model. The model utilizes a planar π-periodic von-Mises distribution for the fiber density distribution ρ(r, M ) around a mean direction M , in the same sense as Gasser et al. [28]. r is the unit orientation vector of a micro-fiber and t := F r is the Eulerian counterpart of the Lagrangian fiber vector. The affine-stretch of a single fiber in the orientation direction r reads The macroscopic free energy corresponding to one family of fibers with the mean direction M and with n fibers per unit volume is defined as where ψ f ib is the free energy function associated with the orientation direction r and | S |= 4π for a unit-sphere. For the isotropic ground matrix, they utilized the neo-Hookean free energy function. The anisotropic free energy function is as follows PSfrag replacements Fig. 12 The polar plots of density distribution (DD) and the directional stiffness (DS) for the HNORS model [37]; (a) in-plane distribution of two families of fibers based on ET dataset for AAA tissue [64], (b) out-of-plane distribution of one families of fibers based on ET dataset for AAA tissue [64] (c) in-plane distribution of one family of fibers based on UT dataset for the linea alba [17], (d) out-of-plane distribution of one family of fibers based on UT dataset for the linea alba [17], (e) in-plane distribution of one family of fibers based on UT dataset for the rectus sheath [55], (f) out-of-plane distribution of one family of fibers based on UT dataset for the rectus sheath [55]. Fig. 13 The polar plots of the density distribution (DD) and the directional stiffness (DS) for the AMDM model [2]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT dataset for the rectus sheath [55].

PSfrag replacements
frag replacements Fig. 14 The polar plots of the density distribution (DD) and the directional stiffness (DS) for the ASMD model [3]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT dataset for the rectus sheath [55].
The contribution of each family of fibers to the macroscopic isochoric Kirchhoff stresses can be expressed as a continuous average, including the orientation distribution function, namely The continuous average in (80) and (82) is approximated by where w i i=1,...,m are the weight factors associated with the discrete orientation directions r i i=1,...,m . The polar plots of the density distribution and the directional stiffness for the AMDM model obtained with the identified parameters for each tissue are depicted in Figure 13. The rotationally symmetric π-periodic, normalized von-Mises distribution reads In this regard, the AMDM model can be considered as the AI counterpart of the GOH model.

AI-based Alastrué -Saéz-Martinez-Doblaré (ASMD) model
As an extension of their previous model, Alastrué et al. [3] included the Bingham distribution in their constitutive model. This distribution function exhibits andipodal symmetry and is expressed as where etr(•) ≡ exp(tr(•)), Z is a diagonal matrix with eigenvalues [κ 1 , κ 2 , κ 3 ], Q is orthogonal orientation matrix such that A = Q · Z · Q T and F 0000 (Z) is defined as where 1 F 1 is a confluent hypergeometric function of the matrix argument. The shape of the distribution is controlled by κ 1 , κ 2 , and κ 3 . In this regard, the ASMD model is a slight modification of the AMDM model that utilizes the Bingham distribution, where the density distribution given in (85) replaces π-periodic von-Mises distribution (84) in the equations (80,82). The polar plots of the density distribution and the directional stiffness for the ASMD model obtained with the identified parameters for each tissue are depicted in Figure 14.
PSfrag replacements Fig. 15 The polar plots of the density distribution (DD) and the directional stiffness (DS) for the DBB model [20]; (a) two families of fibers distribution based on ET dataset for AAA tissue [64], (b) one family of fibers distribution based on UT dataset for the linea alba [17], (c) one family of fibers distribution based on UT dataset for the rectus sheath [55].

AI-based Driessen-Bouten-Baaijens (DBB) model
Driessen et al. [20] presented an extended version of the Holzapfel et al. [36] (HGO) model which included a fiber volume fraction. They applied the rule of mixtures and expressed the isochoric Kirchhoff stress for multiple fiber directions as follows where τ m is the isotropic matrix stress, v f is the volume fraction of fibers and τ f is the fiber stress for a given orientation direction r. The isotropic matrix material is modeled as a neo-Hookean material with a shear modulus µ. The stress expressions for the isotropic matrix and fibers are where the fiber stretch λ f is given in (79). The DBB model utilizes a planar Gaussian distribution around the mean fiber orientation for the orientation fiber volume content where ϑ is the mean value, σ is the standard deviation and A is the normalization constant which is defined as Therein, v tot is the total fiber volume fraction. Unlike the AMDM and ASMD models, the DBB model excludes the volume fraction of the matrix for each orientation direction leading to a structure 1 where, the directional stress additively decomposes into the matrix contribution and collagen fibers contribution proportional to their volume fraction in each orientation direction, respectively. The polar plots of the density distribution for the DBB model obtained with the identified parameters for each tissue are depicted in Figure 15.

Experimental Studies
In this section, we outline the experiments that are used for the comparative investigations made in this paper. To this end, three distinct human tissues; (i) the aneurysmatic abdominal aorta, (ii) the linea alba, and (iii) the anterior rectus sheath are selected for the subsequent analysis. Histologically, the AAA tissue has two alternating families of collagen fibers each orientated at an angle approximately 30 • to the circumferential direction [64], the linea alba has collagen fibers situated in the transverse direction [17], while the rectus sheath's principal fiber direction is in the longitudinal direction [55]. In the experimental studies, the aneursymatic arterial tissue is subjected to equibiaxial tension, whereas the data of the linea alba and anterior rectus sheath originate from uniaxial tension tests executed in two mutually orthogonal directions. In this part, the experimental protocol is briefly reviewed.

Uniaxial tension
Uniaxial tensile tests are performed with samples having a length-to-width ratio of at least 4:1, see Figure 4(a). To experimentally investigate anisotropy, samples are usually cut and tested along both the e 1 -and e 2 -directions, as seen in Figure 4, with the assumption that no anisotropy derives from the e 3 -direction. As the reference axis can change depending on the samples being tested, all following discussion will take place in the context of the e 1 -direction samples, with the understanding that the same principles can be easily adapted by the reader for the e 2 -direction samples. Typically for the uniaxial tension tests of metals or polymers, dogbone samples are employed to reduce the influence of end effects, and encourage fracture within the middle of the specimen and not at the grip location. However, for soft biological tissues, rectangular samples are often used instead due to the soft and delicate nature of the material, rendering the punching of precise dogbone samples generally very difficult.
The first study selected to present uniaxial tension data on soft biological tissues was by Cooney et al. [17], who investigated the uniaxial and equibiaxial tensile behavior of the human linea alba. Cooney et al. [17] obtained 13 freshly frozen human cadaveric abdominal walls for their study and from these extracted the linea alba for their experiments. The linea alba is a collagenous part of the ventral abdominal wall whose collagen fibers are known to be anisotropically arranged. Prior to sample extraction, the human cadaveric abdominal walls were allowed to defrost for 36 hours at 4 • C. The linea alba were then extracted from each abdominal wall and stored in a phosphatebuffered saline solution until testing. For the uniaxial tension tests, 14 rectangular samples in total were obtained from 7 of the linea alba specimens. Seven samples were cut in both the longitudinal and transverse directions, with longitudinal and transverse referring to their respective anatomical axes. The length-to-width ratio of the rectangular samples was as close to the uniaxial test condition as possible for the size of available tissue, with this equaling approximately 2:1 (length:width) for both directions. Prior to testing, samples were mounted in grips lined with emery paper to reduce any slippage during loading. The grips were tightened using a 0.2 N-m torque wrench to prevent over-tightening of the grips causing damage to the tissue. Six black dots were applied to the surface of the samples, before securing them in the modular mechanical testing machine which was fitted with a 300 N load cell. During testing, a high definition camera was used to collect images at a rate of 2 Hz. The post-processing of these images allowed for the calculation of the sample strains by analyzing the deformation of each black dot throughout the stretching. A prestress of 0.1 MPa was applied to all samples to remove any slack present. The tests were conducted at a quasistatic strain rate of 28.5% min −1 . The transverse direction was found to be much stiffer than the longitudinal direction, displaying the anisotropy of the tissue. Cooney et al. [17] found the slope of the most linear section of the curve to be approximately 72 MPa in the transverse direction and 8 MPa in the longitudinal direction. The second study selected was by Martins et al. [55], who investigated the uniaxial tensile behavior of the human anterior rectus sheath. Martins et al. extracted the tissue samples from 12 fresh female cadavers. Six rectangular samples were cut both parallel (longitudinal) and perpendicular (transverse) to the fiber direction with a length-to-width ratio of 4:1. The thickness of the samples were 1.00 ± 0.17 mm. A combination of sandpaper and Velcro tape were used between the grips to stop the soft tissue sample from slipping during testing. A 200 N load cell was fitted to the traction machine to obtain the force, and the displacement of the samples was measured using a displacement sensor. The uniaxial tests were carried out at room temperature and were conducted until rupture at a rate of 5 mm min −1 . The tensile behavior was observed to be nonlinear in both the fiber direction and transverse to the fiber direction.

Equibiaxial tension
Equibiaxial tensile tests are utilized in the characterization of soft biological tissues as the test condition is often considered to be more representative of the stresses experienced in vivo [17]. Square samples are usually preferred for the tests and the samples are equally loaded under tension on all four sides, as seen in Figure 4(c). Due to this, the principal stretch in the e 1 -and e 2 -directions are equal to the experimental stretch, i.e. λ 1 = λ 2 = λ, leaving the e 3 -direction the only one unconstrained. If the material is considered incompressible, and due to the assumption of symmetry, λ 3 = λ −2 . The equibiaxial tension study selected was an investigation by Niestrawska et al. [64] on human abdominal aortas. They conducted equibiaxial tension experiments on both healthy and aneurysmatic abdominal aortas. The data used here was extracted from the experimental results on the abdominal aortic aneurysms which were obtained through surgical procedures conducted at the Department of Vascular Surgery, Medical University Graz, Austria. Aneurysms are defined as irreversible, localized dilatation of a vessel, which can in some cases lead to complete wall rupture. Eleven wall samples in total were collected from open aneurysm repair at the anterior side. The samples were stored at 4 • C in Dublecco's modified Eagle's medium prior to testing. To prepare for testing, 20 × 20mm patches were cut from the aneurysmatic abdominal aorta walls. Some were large enough for two test specimens to be prepared from a single extracted sample. The specimens were tested with their layers intact as the authors found a clear separation of layers of the aneurysms to be impossible. At this point, the thickness of the specimens were measured. A scatter pattern of black dots was then applied to the surface of the specimens using a spray to allow the displacements to be optically measured. To mount the specimens in the biaxial testing machine, hooked surgical sutures were used. The specimens were then submerged in a 0.9% physiological saline bath which was heated up to 37 ± 0.1 • C. The tests were conducted using a stretch-driven protocol, starting at 2.5% deformation and increasing in 0.025 stretch steps until rupture. This was carried out at a quasi-static strain rate of 3mm min −1 , and with a stretch ratio of λ axial : λ circ = 1 : 1, where λ axial is the stretch in the axial direction and λ circ is the stretch in the circumferential direction. Four preconditioning cycles were conducted after each increase in step, with the fifth cycle being recorded for data analysis. The results presented in the study by Niestrawska et al. [64] included the Cauchy stress-stretch data in both the axial and circumferential directions of 12 patch specimens. The results displayed large variability in mechanical response. However, despite the variability, when comparing the Cauchy stress at 1.15 stretch, the stresses in the circumferential direction were consistently higher than the axial direction. The experimental data chosen here was from the specimen with the behavior that was approximately median out of all the specimens (specimen AAA-1.2).

Parameter identification procedure
The parameter optimization procedure introduced by Dal et al. [18] has been adopted for the constitutive models selected here for anisotropic soft tissues. The parameter identification process is conducted based on error expressions for the uniaxial tension experiments in the e 1 -direction, uniaxial tension experiments in the e 2 -direction, and the equibiaxial tension experiment in the e 1 − e 2 directions, respectively, where P 11 and P 22 are the first Piola-Kirchhoff stresses, and n UT1 , n UT2 , n ET are number of data points for the UT in e 1 , UT in e 2 , and ET experiments, respectively. The total cost function for the UT tests and the ET test are The total cost functions are presented for the UT tests and ET tests individually since there is a lack of UT and ET test data that belong exactly to the same tissue. The parameters are extended to include the weights in 92, ξ := {ζ, w} which is obtained from the minimization principle along with the admissible parameter domain where w 1 + w 2 = 1. The domain D is the physically admissible domain for the material parameters ζ. The gradient type optimization is conducted by the Fmincon function in Matlab where the extended cost function is applied to equality constraints h i (ξ) and inequality constraints g i (ξ), respectively. To obtain an optimum solution, the variation of cost functions with respect to the extended parameters requires inclusive of the Karush-Kuhn-Tucker optimality conditions for inequality constraints where λ ine i are the Lagrange multipliers for the inequality constraint. The penalty parameters λ eq i enforce the equality constraint, The Fmincon function in Matlab is used to minimize E, subject to the linear equality A eq ξ = b eq and inequality Aξ ≤ b. Therein, A eq is the coefficient matrix for the equality constraint, b eq is the vector for the equality constraint, A is the coefficient matrix for the inequality constraints, b is the vector for the inequality constraint, and ξ 0 are the initial points.
In this study, a hybrid optimization procedure is employed in order to reach the best parameter set. We utilize the above outlined minimization principle in conjunction with the genetic algorithm presented in Dal et al. [18]. The material parameter space is first trained with the genetic algorithm and the best parameter sets resulting from the genetic algorithm are used as starting points for the gradient search algorithm.

Comparison of hyperelastic anisotropic models
The quality of fit metric (χ 2 ) is used to compare the fitting performance of models. The quality of fit parameter for the uniaxial dataset of Cooney et al. [17], and Martins et al. [55] is where n UT1 and n UT2 are the number of data points and P UT1

11
and P UT2

22
are the first Piola-Kirchhoff stresses for the UT test in the e 1 and e 2 directions, respectively. Similarly, for the equibiaxial loading case, the quality of fit parameter for the fitting of the Niestrawska et al. [64] dataset is where n ET1 and n ET2 are the number of data points and P ET1

11
and P ET2

22
are the first Piola-Kirchhoff stresses for the ET test in the e 1 − e 2 directions. The quality of fit metric for each model has been presented in three regions based on the stretch ranges, g 6 Results and discussion Table 1 Models sorted according to the quality of fit to the equibiaxial dataset of AAA tissue [64], uniaxial dataset of the linea alba [17] and uniaxial dataset of the rectus sheath [55]. The models are compared based on their quality of fit metric using the ET dataset for AAA tissue [64], UT dataset for the linea alba [17], and UT dataset for the rectus sheath [55].  represent the simultaneous fit results of the models to each dataset. In the parameter optimization procedure, we have adhered to the histological information provided by Niestrawska et al. [64], Cooney et al. [17], and Martins et al. [55] in that, for each model, we have imposed the same mean fiber directions. However, because of the lack of data on the fiber dispersion, we decided to treat the distribution parameters as model parameters to be determined during the parameter identification procedure. In the ideal case, the distribution parameters should be obtained through fitting with histologically obtained fiber density distribution data. For nine anisotropic models, identified parameters and error bounds for AAA tissue, linea alba, and rectus sheath can be found in Appendix.
AAA tissue: Using the ET dataset for AAA tissue [64], the models are sorted regarding the quality of fit metric, and the results are listed in Table 1(column #1). The ET dataset for AAA tissue was obtained from an arterial wall specimen with two families of fibers. All models except the HGO [36], and the OS [65] models fit the AAA tissue quite successfully. The OS model fails to capture the transition of the stress-stretch relation from a linear phase to an exponential phase for the curves in both the axial and circumferential directions. On the other hand, the HGO model fail to simultaneously predict the equibiaxial stress-stretch behavior in both the circumferential and axial directions whereas they are capable of fitting either the circumferential or axial stressstretch behavior successfully. The first six models showed a remarkable fitting performance to the ET dataset of AAA tissue [2,3,20,28,37,41,58]. It can be easily seen that the models taking the fiber dispersion into account are relatively more successful. For instance, HSGR and NY models use a phenomenological parameter to represent a measure of dispersion. The most successful model according to the fitting performance of the AAA tissue ET dataset is the six-parameter HNORS model [37] based on the bivariate von-Mises distribution that considers the in-and out-of-plane dispersion of fibers. The five-parameter HSGR model, and the four-parameter NY model [58] also exhibit a remarkable fitting performance with fewer material parameters. AI-based AMDM model [2] and GST-based GOH model [28] also demonstrated excellent fitting performance; these two models consider fiber dispersion with an identical von-Mises distribution function and use an equivalent free energy function. The only difference between these models is the integration approach. Their fitting performance is almost equivalent, however, the computational cost of the GOH model [28] is significantly less compared to the AMDM model [2]. The computational cost of the AI-based formulations are considerably higher than the GST-based ones, see also [39].
Linea alba: The comparison result of the models based on the UT dataset for the linea alba is listed in Table 1(column #2). The UT dataset for linea alba exhibits a weak exponential form of the stress-stretch behavior. In line with the comparison for the AAA tissue dataset, fiber dispersion models have a better fitting performance for the linea alba dataset compared to invariant-based models. This is particularly due to the dispersion of fibers of the human linea alba for which the collagen fibers have a less sophisticated distribution than the collagen distribution within the arterial wall. Thus, using higher-order distribution functions to represent fiber distributions of the linea alba loses importance and increases computational complexity. For the linea alba, the HNORS model [38] exhibits the best fitting performance.
Rectus sheath: The rectus sheath tissue exhibits a relatively softer mechanical response with a weaker degree of anisotropy. The models are sorted according to the quality of fit metric and are listed in Table 1(column #3). It can be observed that in both directions, the stretch-stress response of the rectus sheath tissue is nonlinear. The NY model [58] has the best fitting performance for the rectus sheath dataset [55]. When considering all 9 models, the NY model [58] is the only model which does not decompose the free energy function into the isotropic part and anisotropic part. Thus, even when the experimental data presents stretch-stress information in the transverse direction to the fibers, the NY model [58] considers this data as fibrous. Therefore, although the NY model [58] does not represent the histological structure accurately, it is able to capture the stretch-stress response mathematically. However, the disadvantage of not decomposing the free energy function can be seen when considering the fitting performance of the NY model [58] on the UT dataset for the linea alba. For the linea alba [17], the data presented in transverse to fiber direction is almost linear, and the NY model [58] cannot fit this data. The next best 5 models for the rectus sheath dataset are fiber dispersion-based models. This could be due to the highly dispersed fiber structure of the rectus sheath; even transverse to the fiber direction, the experimental data shows a J-shape response. Transversely isotropic models such as the HGO model cannot capture this mechanical response, since the non-fibrous part of the tissue is modeled by the neo-Hookean model. However, another transversely isotropic model, the OS model [65] is able to show better fitting performance than the HGO model. Similar to the case of the NY model, the performance of the OS model is not due to the accurate representation of the histological structure of the tissue. The reason is that OS model uses the Gent model [30] for the non-fibrous part of a tissue, which is a non-linear model. Furthermore, GST-based dispersion models have weak performance on rectus sheath data. GST-based dispersion models such as GOH [28] and HNORS [37] models impose tension only condition for mean fiber directions. However, in the case of rectus sheath, the mean direction of fibers is along longitudinal direction and as a result GOH [28] and HNORS [37] models are reduced to neo-Hookean model to predict mechanics behavior along transverse direction to the fibers. The same behavior can also be seen in HGO [36], HSGR [41] models.
Summary: The comparison of the results demonstrate that the fiber dispersion-based models have a superior fitting performance over the strain invariant-based models. GST-based dispersion models are more cost effective compared to the AI-based dispersion formulations regarding the computational time. On the other hand, the quality of fit of AI-and GST-based formulations are comparable for AAA tissue, see for example the AMDM [2] and GOH [28] models. However, GSTbased approaches have poor performance on tissues with single family of fibers due to tension only condition. Moreover, the fitting performance of the AI-based formulations for the three distinct tissues shows that these formulations are able to adapt themselves for a variety of tissues having different families of fiber architecture and distribution. The OS model [65], based on the functional form of the Gent model originally proposed for the non-Gaussian chain statistics of the rubber network, is not as successful as the models based on the Fung-type exponential representations for the anisotropic part of the free energy function, see for example HSGR [41], GOH [28], and AMDM [2] models.