# Integrity basis of polyconvex invariants for modeling hyperelastic orthotropic materials - Application to the mechanical response of passive ventricular myocardium

Renye Cai, Frédéric Holweck, Zhi-Qiang Feng, François Peyraut

► To cite this version:

Renye Cai, Frédéric Holweck, Zhi-Qiang Feng, François Peyraut. Integrity basis of polyconvex invariants for modeling hyperelastic orthotropic materials - Application to the mechanical response of passive ventricular myocardium. International Journal of Non-Linear Mechanics, 2021, 133, pp.103713. 10.1016/j.ijnonlinmec.2021.103713. hal-03204826

HAL Id: hal-03204826

https://hal.science/hal-03204826

Submitted on 5 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Integrity basis of polyconvex invariants for modeling hyperelastic orthotropic materials — Application to the mechanical response of passive ventricular myocardium

Renye Cai<sup>a</sup>, Frédéric Holweck<sup>b</sup>, Zhi-Qiang Feng<sup>c,d</sup>, François Peyraut<sup>b,*</sup>

The present paper proposes a new Strain Energy Function (SEF) for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. In order to build our SEF, we have followed a classical strategy based on exponential functions, but we have chosen to work with polyconvex invariants instead of the standard ones. Actually, in the context of hyperelastic problems, the polyconvexity of the strain energy density is considered as a prerequisite for ensuring the existence of solutions. By selecting a set of polyconvex invariants, we demonstrate that our model can predict the experimental data with 6 different shear modes applied to passive ventricular myocardium.

## 1. Introduction

Understanding the behavior of anisotropic hyperelastic materials is of major importance for scientists because their modeling has a wide range of applications in engineering biosciences such as in health therapeutic, medical prosthesis, ergonomics or virtual surgery. The mechanical study of the shear deformation of myocardial layers is for example useful because these deformations are considered to play an important role in the mechanical behavior of the heart [1].

These past ten years, many works have been performed to investigate the structurally based model originally proposed by Holzapfel and Ogden in [2] in relation to the tests carried out in [1] where the orthotropic nature of the ventricular myocardium has been proven. Erikson et al. [3] have for example improved this model with structure parameters allowing the quantification of the degree of dispersion based on measured fiber and sheet angle data. Melnik et al. [4] have modified the hyperelastic Holzapfel-Ogden model by using Generalized Structure Tensors (GSTs) in order to account for fiber dispersion in every term of the strain-energy function. McEvoy et al. [5] have studied the effect of compressibility by combining a nonlinear volumetric hyperelastic component with the Holzapfel-Ogden anisotropic hyperelastic component for myocardium fibers. Palit et al. [6] have implemented the Holzapfel-Ogden constitutive law in the MSC-Marc finite element software in order to inversely estimate the constitutive parameters of the model. Very recently, Li et al. [7] have extended the

Holzapfel-Ogden model by accounting for the mixed-invariants in the fiber-normal and sheet-normal directions.

All these models, like in most of the papers published in the literature, separate the SEF into an isotropic part and an anisotropic part. The first part is used to model the low strain behavior of the ground matrix and the second part accounts for the behavior of the fibers at higher strain. A first alternative to this standard approach is to introduce free-invariants SEF [8]. Another alternative consists in using an integrity basis of invariants as proposed in [9] and mixing them in a single SEF [10]. This second alternative, mathematically justified by the theory of invariant polynomials and by Noether's theorem, is inspired by the pioneering work of Thionnet et al. [11]. One advantage of this second alternative is to work with polyconvex invariants, the polyconvexity of the strain energy density being considered as a prerequisite for ensuring the existence of solutions in compatibility with physical requirements [12]. A wide survey with many proofs on polyconvexity of isotropic and transversely isotropic functions can be viewed in [13].

The present paper proposes a new SEF for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. In order to build our SEF, we have followed a strategy based on exponential functions, as proposed by Holzapfel et al. [2], but we have selected polyconvex invariants instead of the standard ones generally used in the literature [2-7]. Working with the set of polyconvex invariants exhibited in [10] allows to replace the classical mixed invariant $I_8$, which is

<footnote><sup>a</sup> School of Automobile and Transportation Engineering, Guangdong Polytechnic Normal University, Guangdong, China</footnote>

<footnote><sup>b</sup> ICB, UMR 6303, CNRS, Univ. Bourgogne Franche-Comté, UTBM, 90010 Belfort, France</footnote>

<footnote><sup>c</sup> School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, China</footnote>

<footnote><sup>d</sup> LMEE Univ-Evry, Université Paris-Saclay, 91020 Evry, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: francois.peyraut@utbm.fr (F. Peyraut).</footnote>
non polyconvex (the proof is in Section 4), by the polyconvex invariant $L_4$ defined by Eq. (10). This replacement provides consistent numerical results with [2] and our model can perfectly match the experimental data obtained by Dokos et al. [1] with 6 different shear modes applied to passive ventricular myocardium. In addition, sufficient conditions ensuring the convexity of the SEF are given in Section 7. Finally, our model can represent stretch effect in the 3 specific directions of passive ventricular myocardium, the myocyte axis direction, the direction lying within the muscle layer and transverse to the myocyte axis, and the direction normal to the muscle layer.

## Notations

A bold-face Latin lowercase letter, say $\mathbf{a}$, and a bold-face Latin capital letter, say $\mathbf{A}$, will denote a vector and second-order tensor, respectively. The standard Euclidean inner product $\langle \cdot, \cdot \rangle$ in a $n$ vector space dimension, its related norm $\| \cdot \|$, and the product $\otimes$. between two vectors $\mathbf{a}$ and $\mathbf{b}$, are respectively defined by:

$$ \langle \mathbf{A}\mathbf{a}, \mathbf{a} \rangle = \sum_{i=1}^{n} A_{ij} a_i a_j; \quad \| \mathbf{a} \| = \sqrt{\sum_{i=1}^{n} (a_i)^2}; \quad (\mathbf{a} \otimes \mathbf{b})_{ij} = a_i b_j $$

## 2. Kinematics and basic continuum mechanics

Consider a continuum body $\mathbf{V}$ with particle $P \in \mathbf{V}$ which is embedded in the three-dimensional Euclidean space at a given instant $t$ of time. As the continuum body $\mathbf{V}$ moves in space from one instant of time to another it occupies a continuous sequence of geometrical regions denoted by $\Omega_0...\Omega_i, \Omega_0$ is referred to as the fixed reference (undeformed or Lagrangian) configuration ($t = 0$) of the body $\mathbf{V}$ while the configuration $t$ is called the current (deformed or Eulerian) configuration. The position $\mathbf{x}$ of particle $P$ in the current (or deformed) configuration can be deduced by the motion $\varphi$ that depends on its initial position vector $\mathbf{X}$ and the time $t$:

$$ \mathbf{x} = \varphi(\mathbf{X}, t) \qquad (1) $$

The deformation gradient matrix $\mathbf{F}$ is defined by:

$$ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \mathbf{I} + \frac{\partial U(\mathbf{X}, t)}{\partial \mathbf{X}} \qquad (2) $$

$\mathbf{I}$ is the unity tensor and $\mathbf{U}(\mathbf{X}, t) = \mathbf{x} - \mathbf{X}$ is the displacement vector of the particle. The tensors $\mathbf{C}$ and $\mathbf{B}$ are the so-called right and left Cauchy-Green strain tensors:

$$ \mathbf{C} = \mathbf{F}^T \mathbf{F}; \quad \mathbf{B} = \mathbf{F} \mathbf{F}^T \qquad (3) $$

In terms of stress tensor, we remind that the second Piola-Kirchhoff stress tensor $\mathbf{S}$ and the corresponding Cauchy stress tensor $\boldsymbol{\sigma}$ are obtained by differentiating a SEF $\mathbf{W}$ with respect to $\mathbf{C}$:

$$ \mathbf{S} = 2 \frac{\partial \mathbf{W}}{\partial \mathbf{C}} - p \mathbf{C}^{-1}; \quad \boldsymbol{\sigma} = \mathbf{J}^{-1} \mathbf{F} \mathbf{S} \mathbf{F}^T \qquad (4) $$

where the extra pressure $p$ is included in the formulation to account for the incompressibility condition $\mathbf{J} = \det(\mathbf{F}) = 1$. Combining the two equalities of Eq. (4) yields to:

$$ \boldsymbol{\sigma} = 2 \mathbf{J}^{-1} \mathbf{F} \frac{\partial \mathbf{W}}{\partial \mathbf{C}} \mathbf{F}^T - p \mathbf{I} \qquad (5) $$

Note that, in the following, the time $t$ will be omitted because we are only interested in static problems.

## 3. Integrity basis of polyconvex invariants for modeling orthotropic materials

The mechanical response of passive ventricular myocardium to simple shear loading is sensitive to the shear mode as described in [1]. In this context, the fiber direction $\mathbf{f}_0$, the direction $\mathbf{s}_0$ transverse to the fiber direction in plane of sheets, and the normal $\mathbf{n}_0$ to sheets play a key role to define the material symmetry group (Fig. 1). In particular, the 3 planes perpendicular to the three directions $\mathbf{f}_0, \mathbf{s}_0$ and $\mathbf{n}_0$ are planes

of symmetry for the material. The material properties remain therefore invariant under the action of 3 reflections $\mathbf{R}(\mathbf{f}_0), \mathbf{R}(\mathbf{s}_0)$ and $\mathbf{R}(\mathbf{n}_0)$ related to these 3 planes; 3 rotations $\mathbf{Q}_\pi(\mathbf{f}_0), \mathbf{Q}_\pi(\mathbf{s}_0)$ and $\mathbf{Q}_\pi(\mathbf{n}_0)$ by an angle $\pi$ around $\mathbf{f}_0, \mathbf{s}_0$ and $\mathbf{n}_0$, and under $\mathbf{I}$ and $-\mathbf{I}$:

$$ \mathbf{R}(\mathbf{f}_0) = \mathbf{I} - 2\mathbf{f}_0 \otimes \mathbf{f}_0 = -\mathbf{Q}_\pi(\mathbf{f}_0) : \mathbf{R}(\mathbf{s}_0) = \mathbf{I} - 2\mathbf{s}_0 \otimes \mathbf{s}_0 = -\mathbf{Q}_\pi(\mathbf{s}_0) \qquad (6) $$

$$ \mathbf{R}(\mathbf{n}_0) = \mathbf{I} - 2\mathbf{n}_0 \otimes \mathbf{n}_0 = -\mathbf{Q}_\pi(\mathbf{n}_0) \qquad (7) $$

The set $S_8 = \{\mathbf{R}(\mathbf{f}_0), \mathbf{R}(\mathbf{s}_0), \mathbf{R}(\mathbf{n}_0), \mathbf{Q}_\pi(\mathbf{f}_0), \mathbf{Q}_\pi(\mathbf{s}_0), \mathbf{Q}_\pi(\mathbf{n}_0), \mathbf{I}, -\mathbf{I}\}$ built with the 8 invariant matrix operators defined by Eqs. (6)–(7) is stable under multiplication, contains the neutral element $\mathbf{I}$ and each element has an inverse in $S_8$. This confers a group structure to $S_8$. Using a mathematical argument based on the Reynolds operator and on the Noether's theorem, Ta et al. [9] have presented a systematic method to find a set of invariants associated with $S_8$ and demonstrated that the following 7 polynomial invariants $K_i$ form an integrity basis of the ring of invariant polynomials under the material symmetry group $S_8$:

$$ \begin{aligned} K_1 &= \rho_1; & K_2 &= \rho_2; & K_3 &= \rho_3; & K_4 &= \rho_4^2; & K_5 &= \rho_5^2; & K_6 &= \rho_6^2; & K_7 &= \rho_4 \rho_5 \rho_6 \end{aligned} \qquad (8) $$

$$ \begin{cases} \rho_1 = \langle \mathbf{C} \mathbf{f}_0, \mathbf{f}_0 \rangle; & \rho_2 = \langle \mathbf{C} \mathbf{s}_0, \mathbf{s}_0 \rangle; & \rho_3 = \langle \mathbf{C} \mathbf{n}_0, \mathbf{n}_0 \rangle \\ \rho_4 = \langle \mathbf{C} \mathbf{f}_0, \mathbf{s}_0 \rangle; & \rho_5 = \langle \mathbf{C} \mathbf{f}_0, \mathbf{n}_0 \rangle; & \rho_6 = \langle \mathbf{C} \mathbf{s}_0, \mathbf{n}_0 \rangle \end{cases} \qquad (9) $$

Later, Cai et al. [10] combined the invariants $K_i$ in order to build an integrity basis of 7 polyconvex invariants $L_i$:

$$ \left\{ \begin{aligned} L_1 &= K_1; & L_2 &= K_2; & L_3 &= K_3; & L_4 &= (K_1 + K_2)^2 + 4K_4; \\ L_5 &= (K_1 + K_3)^2 + 4K_5 \\ L_6 &= (K_2 + K_3)^2 + 4K_6; \\ L_7 &= K_7 + \frac{1}{2} (K_1 K_2 K_3 - K_1 K_6 - K_2 K_5 - K_3 K_4) \end{aligned} \right. \qquad (10) $$

It is also demonstrated in [10] that $L_7$ contributes in Eq. (5) under the form of a redundant term with the extra pressure $p$. $L_7$ can therefore be removed from the set of invariants. Consequently, in order to benefit from the property of polyconvexity, we assume that the SEF $\mathbf{W}$ depends on the six first polyconvex invariants $L_i$. Therefore, Eq. (5) yields to:

$$ \boldsymbol{\sigma} = 2 \mathbf{J}^{-1} \mathbf{F} \left( \sum_{i=1}^{6} \frac{\partial \mathbf{W}}{\partial L_i} \frac{\partial L_i}{\partial \mathbf{C}} \right) \mathbf{F}^T - p \mathbf{I} \qquad (11) $$

The derivatives $\frac{\partial L_i}{\partial \mathbf{C}}$ embedded in Eq. (11) are calculated straightforwardly from Eqs. (8)–(10):

$$ \frac{\partial L_1}{\partial \mathbf{C}} = \mathbf{f}_0 \otimes \mathbf{f}_0; \quad \frac{\partial L_2}{\partial \mathbf{C}} = \mathbf{s}_0 \otimes \mathbf{s}_0; \quad \frac{\partial L_3}{\partial \mathbf{C}} = \mathbf{n}_0 \otimes \mathbf{n}_0 \qquad (12) $$

$$ \frac{\partial L_4}{\partial \mathbf{C}} = 2 \left\{ (L_1 + L_2) (\mathbf{f}_0 \otimes \mathbf{f}_0 + \mathbf{s}_0 \otimes \mathbf{s}_0) + \sqrt{L_4 - (L_1 + L_2)^2} (\mathbf{f}_0 \otimes \mathbf{s}_0 + \mathbf{s}_0 \otimes \mathbf{f}_0) \right\} \qquad (13) $$

$$ \frac{\partial L_5}{\partial \mathbf{C}} = 2 \left\{ (L_1 + L_3) (\mathbf{f}_0 \otimes \mathbf{f}_0 + \mathbf{n}_0 \otimes \mathbf{n}_0) + \sqrt{L_5 - (L_1 + L_3)^2} (\mathbf{f}_0 \otimes \mathbf{n}_0 + \mathbf{n}_0 \otimes \mathbf{f}_0) \right\} \qquad (14) $$

$$ \frac{\partial L_6}{\partial \mathbf{C}} = 2 \left\{ (L_2 + L_3) (\mathbf{s}_0 \otimes \mathbf{s}_0 + \mathbf{n}_0 \otimes \mathbf{n}_0) + \sqrt{L_6 - (L_2 + L_3)^2} (\mathbf{s}_0 \otimes \mathbf{n}_0 + \mathbf{n}_0 \otimes \mathbf{s}_0) \right\} \qquad (15) $$

In order to calculate the derivatives $\frac{\partial \mathbf{W}}{\partial L_i}$ also included in Eq. (11), it is necessary to define how the SEF $\mathbf{W}$ depends on the invariants $L_i$. This issue is discussed in the next section.

## 4. Strain energy function

Dokos et al. [1] have demonstrated that passive ventricular myocardium behaves differently depending on the plane of shear. The
<figure><img src="image_2.png" /><figcaption>Fig. 1. 6 possible modes of simple shear test in different planes.</figcaption></figure>

6 different modes of shear are shown in Fig. 1 where the myocyte axis representing the fiber direction is depicted with green full lines. Building a consistent model thus requires accounting for the myocyte axis direction $f_0$, the direction $s_0$ lying within the muscle layer and transverse to $f_0$, and the direction $n_0$ normal to the muscle layer.

Due to the exponential trends of the stress observed experimentally, Holzapfel et al. [2] have proposed a SEF made of exponential functions:

$$ W = \frac{a}{2b} \exp[b(I_1 - 3)] + \sum_{i=f,s} \frac{a_i}{2b_i} \{ \exp[b_i(I_{4i} - 1)^2] - 1 \} + \frac{a_{fs}}{2b_{fs}} \{ \exp[b_{fs}I_{8fs}^2] - 1 \} \quad (16) $$

where $a, b, a_f, a_s, b_f, b_s, a_{fs}$ and $b_{fs}$ are 8 positive material constants, $I_1$ is the first principal invariant of $C$, and $I_{4i}$, $I_{4s}$ and $I_{8fs}$ are the classical mixed invariants combining $C$ with $f_0$ and $s_0$:

$$ I_1 = \mathrm{Tr}(C); \quad I_{4f} = \langle Cf_0, f_0 \rangle; \quad I_{4s} = \langle Cs_0, s_0 \rangle; \quad I_{8fs} = \langle Cf_0, s_0 \rangle \quad (17) $$

It is easy to check that the standard invariants (17) are connected to the polyconvex invariants $L_i$ introduced by Eqs. (8)–(10) by:

$$ I_1 = L_1 + L_2 + L_3; \quad I_{4f} = L_1; \quad I_{4s} = L_2; \quad I_{8fs} = \frac{1}{2}\sqrt{L_4 - (L_1 + L_2)^2} \quad (18) $$

And conversely:

$$ L_1 = I_{4f}; \quad L_2 = I_{4s}; \quad L_3 = I_1 - I_{4f} - I_{4s}; \quad L_4 = 4I_{8fs}^2 + (I_{4f} + I_{4s})^2 \quad (19) $$

However, it should be also noted that, if the classical invariants $I_1, I_{4f}$ and $I_{4s}$ are polyconvex, as convex functions of $F$, it is not the case for $I_{8fs}$. Actually, from Eq. (17), the second derivative of $I_{8fs}$ with respect to $F$ is:

$$ \frac{d^2 I_{8fs}}{dF^2} (F)(dF)(dG) = \langle dFf_0, dGs_0 \rangle + \langle dFs_0, dGf_0 \rangle \quad (20) $$

$$ \Rightarrow \frac{d^2 I_{8fs}}{dF^2} (F)(dF)(dF) = 2 \langle dFf_0, dFs_0 \rangle \quad (21) $$

In order to demonstrate that $I_{8fs}$ is non polyconvex, we just have to find a counterexample with a matrix $dF$ giving a negative value in Eq. (21). We choose this $dF$ as follows:

$$ dF = -f_0 \otimes f_0 + f_0 \otimes s_0 + s_0 \otimes f_0 \quad (22) $$

One reports Eq. (22) in the inner product of Eq. (21):

$$ \langle dFf_0, dFs_0 \rangle = \langle -f_0 + s_0, f_0 \rangle = \langle -f_0, f_0 \rangle = -1 < 0 \quad (23) $$

The proof of the non-polyconvexity of the classical mixed invariant $I_{8fs}$ is complete. Besides, it is underlined in [2] that $I_{8fs}$ generally has a destabilizing influence. It turns out that if polyconvexity is a prerequisite for building SEF in order to ensure existence of solutions, it would be better to choose the set of invariants defined by Eqs. (8)–(10) rather than the ones defined by Eq. (17). In order to make further comparison between the model introduced in [2] and our model, it is useful at this stage to remind the Cauchy stress calculated in [2]:

$$ \sigma = a \exp[b(I_1 - 3)] \mathbf{B} - p\mathbf{I} + 2a_f (I_{4f} - 1) \exp[b_f(I_{4f} - 1)^2] f \otimes f + 2a_s (I_{4s} - 1) \exp[b_s(I_{4s} - 1)^2] s \otimes s + a_{fs} I_{8fs} \exp[b_{fs}I_{8fs}^2] (f \otimes s + s \otimes f) \quad (24) $$

where $f$ and $s$ represent the deformed directions, respectively for the fiber and for the direction transverse to the fiber direction in plane of sheets:

$$ f = Ff_0; \quad s = Fs_0; \quad n = Fn_0 \quad (25) $$

To build our model, we proceed at the beginning in the same manner as in [2], by expressing the SEF in the form of exponential functions such as in Eq. (16). However, instead of using the classical invariants introduced by Eq. (17), we select the set of polyconvex invariants defined by Eqs. (8)–(10):

$$ W = \sum_{i=1}^{3} \frac{\alpha_i}{2\rho_i} \{ \exp[\beta_i (L_i - 1)^2] - 1 \} + \frac{\alpha_4}{2\rho_4} \{ \exp[\beta_4 (L_4 - 4)^2] - 1 \} \quad (26) $$

The 8 coefficients $\alpha_1, \alpha_2, \alpha_3, \alpha_4, \beta_1, \beta_2, \beta_3$ and $\beta_4$ are the material parameters of our model, the $\alpha_i$ parameters having the unit of stress while the parameters $\beta_i$ are dimensionless. Note that we have omitted in Eq. (26) the dependence of $W$ with respect to the invariants $L_5$ and $L_6$ in order to keep the same number of invariants as in the model proposed in [2]. It is also remarked that the SEF defined by Eq. (26) is null if the material is at rest. Actually, if the displacement is equal to zero, it comes from
Eqs. (2)–(3) that both tensors F and C are equal to the identity matrix I. Therefore, it follows from Eqs. (8)–(10) that: $L_1 = L_2 = L_3 = 1$; $L_4 = 4$. Consequently, $W = 0$. We now calculate the Cauchy stress corresponding to our model by reporting Eqs. (12)–(15) in Eq. (11) and by assuming the incompressible nature of the material (J = 1):

$$ \sigma = 2F \left\{ \left[ \frac{\partial W}{\partial L_1} + 2 \frac{\partial W}{\partial L_4} (L_1 + L_2) \right] f_0 \otimes f_0 + \left[ \frac{\partial W}{\partial L_2} + 2 \frac{\partial W}{\partial L_4} (L_1 + L_2) \right] s_0 \otimes s_0 + \frac{\partial W}{\partial L_3} n_0 \otimes n_0 + 2 \frac{\partial W}{\partial L_4} \sqrt{L_4 - (L_1 + L_2)^2} (f_0 \otimes s_0 + s_0 \otimes f_0) \right\} F^T - pI \quad (27) $$

We next report the well-known property $n_0 \otimes n_0 = I - f_0 \otimes f_0 - s_0 \otimes s_0$ in Eq. (27) and we use Eq. (25):

$$ \sigma = 2 \frac{\partial W}{\partial L_3} \mathbf{B} - p\mathbf{I} + 2 \left[ \frac{\partial W}{\partial L_1} - \frac{\partial W}{\partial L_3} + 2 \frac{\partial W}{\partial L_4} (L_1 + L_2) \right] f \otimes f + 2 \left[ \frac{\partial W}{\partial L_2} - \frac{\partial W}{\partial L_3} + 2 \frac{\partial W}{\partial L_4} (L_1 + L_2) \right] s \otimes s + 4 \frac{\partial W}{\partial L_4} \sqrt{L_4 - (L_1 + L_2)^2} (f \otimes s + s \otimes f) \quad (28) $$

We finally calculate the derivatives $\frac{\partial W}{\partial L_i}$ from Eq. (26):

$$ \begin{aligned} \frac{\partial W}{\partial L_i} &= \alpha_i (L_i - 1) \exp \left[ \beta_i (L_i - 1)^2 \right], \quad i = 1, 2, 3; \\ \frac{\partial W}{\partial L_4} &= \alpha_4 (L_4 - 4) \exp \left[ \beta_4 (L_4 - 4)^2 \right] \end{aligned} \quad (29) $$

In order to compare the Cauchy stress of Eq. (28) with Eq. (24) extracted from [2], one reports Eq. (29) in Eq. (28):

$$ \sigma = A\mathbf{B} - p\mathbf{I} + [\mathbf{B} - A + C]\mathbf{f} \otimes \mathbf{f} + [D - A + C]\mathbf{s} \otimes \mathbf{s} + E(\mathbf{f} \otimes \mathbf{s} + \mathbf{s} \otimes \mathbf{f}) \quad (30) $$

$$ A = 2\alpha_3 (L_3 - 1) \exp \left[ \beta_3 (L_3 - 1)^2 \right]; \quad B = 2\alpha_1 (L_1 - 1) \exp \left[ \beta_1 (L_1 - 1)^2 \right] \quad (31) $$

$$ C = 4\alpha_4 (L_4 - 4) (L_1 + L_2) \exp \left[ \beta_4 (L_4 - 4)^2 \right]; \quad D = 2\alpha_2 (L_2 - 1) \exp \left[ \beta_2 (L_2 - 1)^2 \right] \quad (32) $$

$$ E = 4\alpha_4 (L_4 - 4) \sqrt{L_4 - (L_1 + L_2)^2} \exp \left[ \beta_4 (L_4 - 4)^2 \right] \quad (33) $$

It is observed that the same 5 tensors B, I, f ⊗ f, s ⊗ s and f ⊗ s + s ⊗ f are involved in both equations (24) and (30). But there exist fundamental differences since Eq. (24) is a sum with a single exponential associated to each term while several exponentials coming from the polyconvex invariants are mixed in the 3 last terms of Eq. (30), according to Eqs. (31)–(33). Beyond working with polyconvex invariants, the interest of Eq. (30) is to account for the invariant $I_{4n} = \langle Cn_0, n_0 \rangle$ and therefore to distinguish explicitly the effect of the stretch in the myocyte axis direction $f_0$, in the direction $s_0$ lying within the muscle layer and transverse to $f_0$, and in the direction $n_0$ normal to the muscle layer. To prove this claim, Eq. (30) is reorganized by factorizing the common terms with A and C:

$$ \sigma = A(B - f \otimes f - s \otimes s) + Bf \otimes f + Ds \otimes s - pI + C(f \otimes f + s \otimes s) + E(f \otimes s + s \otimes f) \quad (34) $$

The first parenthesis in Eq. (34) is simplified by reporting inside the following result which holds in the basis ($f_0$, $s_0$, $n_0$):

$$ \begin{aligned} \mathbf{B} &= f \otimes f + s \otimes s + n \otimes n \\ \Rightarrow \sigma &= An \otimes n + Bf \otimes f + Ds \otimes s - pI + C(f \otimes f + s \otimes s) \end{aligned} \quad (35) $$

$$ + E(f \otimes s + s \otimes f) \quad (36) $$

It is noted from Eqs. (19), (31) and (32) that the coefficients B and D are linked to the invariants $I_{4f}$ and $I_{4s}$:

$$ \begin{aligned} B &= 2\alpha_1 (I_{4f} - 1) \exp \left[ \beta_1 (I_{4f} - 1)^2 \right]; \\ D &= 2\alpha_2 (I_{4s} - 1) \exp \left[ \beta_2 (I_{4s} - 1)^2 \right] \end{aligned} \quad (37) $$

That means that the second and third terms of Eq. (36) represent the effect of the stretch in the directions $f_0$ and $s_0$. Regarding the first term of Eq. (36), it follows from Eqs. (17), (19) and (31) that $L_3$ is equal to the mixed invariant $I_{4n}$:

$$ L_3 = \mathrm{Tr}(C[I - f_0 \otimes f_0 - s_0 \otimes s_0]) = \mathrm{Tr}(Cn_0 \otimes n_0) = I_{4n} \quad (38) $$

The coefficient A is therefore linked to $I_{4n}$ via Eqs. (31) and (38):

$$ A = 2\alpha_3 (I_{4n} - 1) \exp \left[ \beta_3 (I_{4n} - 1)^2 \right] \quad (39) $$

That means that the first term of Eq. (36) represents the effect of the stretch in the direction $n_0$. It is thus possible to give a physical meaning to each of the 3 first terms of Eq. (36) with respect to the specific geometry of the material. This would be useful if one needed to focus specifically on the effect of stretching in the $n_0$ direction. This is a significant difference with the model proposed in [2] which uses the 3 invariants $I_1$, $I_{4f}$ and $I_{4s}$ instead of $I_{4f}$, $I_{4s}$ and $I_{4n}$.

## 5. Simple shear test

A shear experiment is driven by the shear deformation k as shown in Fig. 2. The displacement U, related for example to a simple shear deformation in the ($n_0$, $f_0$) plane, is expressed with respect to k as follows:

$$ U = \begin{Bmatrix} kn_0 \\ 0 \\ 0 \end{Bmatrix} \quad (40) $$

It follows from Eqs. (2), (3) and (40) that:

$$ F = \begin{pmatrix} 1 & 0 & k \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}; \quad C = \begin{pmatrix} 1 & 0 & k \\ 0 & 1 & 0 \\ k & 0 & k^2 + 1 \end{pmatrix}; \quad \mathbf{B} = \begin{pmatrix} k^2 + 1 & 0 & k \\ 0 & 1 & 0 \\ k & 0 & 1 \end{pmatrix} \quad (41) $$

Applying to the 6 different modes of shear the same kind of calculation as the one used for establishing Eqs. (40) and (41), it is easy to obtain the displacement U, the gradient deformation matrix F and the strain tensors C and B related to each of these 6 modes (Table 1).

By using the data contained in Table 1, the expressions of the polyconvex invariants $L_1$, $L_2$, $L_3$ and $L_4$ with respect to k (Table 2) come from Eqs. (8)–(10). Finally, by using Tables 1 and 2, Eq. (28) gives the Cauchy stress expressions corresponding to each of the 6 shear modes (Table 3). Since the (nf) and (ns) modes have the same invariant $L_3$, they also have the same shear stresses. This is consistent with the experimental observations since the responses for the (nf) and (ns) modes are indistinguishable as reported in [1].

## 6. Comparison between the model of Eq. (26) and the experimental data of Dokos et al. [1]

In order to identify the 8 material parameters $\alpha_1, \alpha_2, \alpha_3, \alpha_4, \beta_1, \beta_2, \beta_3$ and $\beta_4$ of the SEF W defined by Eq. (26), we have used the classical coefficient of determination $R^2$:

$$ R^2 = 1 - \frac{S_{\mathrm{res}}}{S_{\mathrm{tot}}} \quad (42) $$

where $S_{\mathrm{res}}$ and $S_{\mathrm{tot}}$ are respectively the residual sum and the total sum of squares over the number of experimental data n:

$$ S_{\mathrm{res}} = \|\mathbf{y} - \mathbf{f}\|^2 = \sum_{i=1}^{n} (y_i - f_i)^2; \quad S_{\mathrm{tot}} = \|\mathbf{y} - \bar{\mathbf{y}}\|^2 = \sum_{i=1}^{n} (y_i - \bar{y})^2 \quad (43) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Simple shear test in the $(n_0, f_0)$ plane.</figcaption></figure>

<table><caption>**Table 1**
Displacement and strain tensors for the 6 shear modes.</caption><thead><tr><th>Shear modes</th><th>U</th><th>F</th><th>C</th><th>B</th></tr></thead><tbody><tr><td>(nf)</td><td><table><tbody><tr><td>kn<sub>0</sub></td><td>(1 0 k)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 0 k)</td></tr><tr><td>0 1 0</td><td>0 1 0</td></tr><tr><td>0 0 1</td><td>k 0 k<sup>2</sup> + 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(k<sup>2</sup> + 1 0 k)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>k</td><td>0 1</td></tr></tbody></table></td></tr><tr><td>(ns)</td><td><table><tbody><tr><td>kn<sub>0</sub></td><td>(1 0 0)</td></tr><tr><td>0</td><td>0 1 k</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 0 0)</td></tr><tr><td>0 1 0</td><td>0 1 k</td></tr><tr><td>0</td><td>k 0 k<sup>2</sup> + 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0</td><td>0 k<sup>2</sup> + 1 k</td></tr><tr><td>0</td><td>k 0 1</td></tr></tbody></table></td></tr><tr><td>(fn)</td><td><table><tbody><tr><td>0</td><td>(1 0 0)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>k f<sub>0</sub></td><td>k 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(k<sup>2</sup> + 1 0 k)</td></tr><tr><td>0 1 0</td><td>0 1 0</td></tr><tr><td>k</td><td>0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 0 k)</td></tr><tr><td>0 1 0</td><td>k k<sup>2</sup> + 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td></tr><tr><td>(fs)</td><td><table><tbody><tr><td>k f<sub>0</sub></td><td>(1 0 0)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>k 1 0</td><td>(k<sup>2</sup> + 1 k 0)</td></tr><tr><td>0 0 1</td><td>0 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(k 0 1)</td></tr><tr><td>k k<sup>2</sup> + 1 0</td><td>0 0 1</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td></tr><tr><td>(sn)</td><td><table><tbody><tr><td>0</td><td>(1 0 0)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>k s<sub>0</sub></td><td>0 k 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 0 0)</td></tr><tr><td>0 k<sup>2</sup> + 1 k</td><td>0 1 0</td></tr><tr><td>0</td><td>k 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 0 0)</td></tr><tr><td>0 1 0</td><td>0 1 k</td></tr><tr><td>k</td><td>0 k<sup>2</sup> + 1</td></tr></tbody></table></td></tr><tr><td>(sf)</td><td><table><tbody><tr><td>k s<sub>0</sub></td><td>(1 k 0)</td></tr><tr><td>0</td><td>0 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>0 1 0</td><td>(1 k 0)</td></tr><tr><td>k k<sup>2</sup> + 1 0</td><td>0 1 0</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td><td><table><tbody><tr><td>k<sup>2</sup> + 1 0 k</td><td>0 1 0</td></tr><tr><td>k 1 0</td><td>0 0 1</td></tr><tr><td>0</td><td>0 0 1</td></tr></tbody></table></td></tr></tbody></table>

<table><caption>**Table 2**
Polyconvex invariants for the 6 shear modes.</caption><thead><tr><th rowspan="2">Shear modes</th><th colspan="4">Polyconvex invariants</th></tr><tr><th>L<sub>1</sub></th><th>L<sub>2</sub></th><th>L<sub>3</sub></th><th>L<sub>4</sub></th></tr></thead><tbody><tr><td>(nf)</td><td>1</td><td>1</td><td>k<sup>2</sup> + 1</td><td>4</td></tr><tr><td>(ns)</td><td>1</td><td>1</td><td>k<sup>2</sup> + 1</td><td>4</td></tr><tr><td>(fn)</td><td>k<sup>2</sup> + 1</td><td>1</td><td>1</td><td>(k<sup>2</sup> + 2)<sup>2</sup></td></tr><tr><td>(fs)</td><td>k<sup>2</sup> + 1</td><td>1</td><td>1</td><td>(k<sup>2</sup> + 2)<sup>2</sup> + 4k<sup>2</sup></td></tr><tr><td>(sn)</td><td>1</td><td>k<sup>2</sup> + 1</td><td>1</td><td>(k<sup>2</sup> + 2)<sup>2</sup></td></tr><tr><td>(sf)</td><td>1</td><td>k<sup>2</sup> + 1</td><td>1</td><td>(k<sup>2</sup> + 2)<sup>2</sup> + 4k<sup>2</sup></td></tr></tbody></table>

<table><caption>**Table 3**
Shear component of the Cauchy stresses for the 6 shear modes.</caption><thead><tr><th>Shear modes</th><th>Cauchy stress - shear component</th></tr></thead><tbody><tr><td>(nf)</td><td>2k <sup>∂W</sup> /∂L<sub>1</sub></td></tr><tr><td>(ns)</td><td>2k <sup>∂W</sup> /∂L<sub>2</sub></td></tr><tr><td>(fn)</td><td>2k [ <sup>∂W</sup> /∂L<sub>1</sub> + 2(k<sup>2</sup> + 2) <sup>∂W</sup> /∂L<sub>2</sub> ]</td></tr><tr><td>(fs)</td><td>2k [ <sup>∂W</sup> /∂L<sub>1</sub> + 2(k<sup>2</sup> + 4) <sup>∂W</sup> /∂L<sub>2</sub> ]</td></tr><tr><td>(sn)</td><td>2k [ <sup>∂W</sup> /∂L<sub>2</sub> + 2(k<sup>2</sup> + 2) <sup>∂W</sup> /∂L<sub>1</sub> ]</td></tr><tr><td>(sf)</td><td>2k [ <sup>∂W</sup> /∂L<sub>2</sub> + 2(k<sup>2</sup> + 4) <sup>∂W</sup> /∂L<sub>1</sub> ]</td></tr></tbody></table>

$y_i$ stands for the experimental data, $f_i$ for the theoretical data and $\bar{y}$ for the mean of the experimental data:

$$ \bar{y} = \frac{1}{n} \sum_{i=1}^{n} y_i \qquad (44) $$

<table><caption>**Table 4**
Identified material parameters.</caption><thead><tr><th>α<sub>1</sub> (kPa)</th><th>α<sub>2</sub> (kPa)</th><th>α<sub>3</sub> (kPa)</th><th>α<sub>4</sub> (kPa)</th></tr></thead><tbody><tr><td>18.877</td><td>2.495</td><td>3.184</td><td>0.168</td></tr><tr><td>β<sub>1</sub>(-)</td><td>β<sub>2</sub>(-)</td><td>β<sub>3</sub>(-)</td><td>β<sub>4</sub>(-)</td></tr><tr><td>19.39</td><td>20.113</td><td>11.543</td><td>0.107</td></tr></tbody></table>

<table><caption>**Table 5**
Coefficient of determination R<sup>2</sup>.</caption><thead><tr><th>Shear modes</th><th>(nf)=(ns)</th><th>(fn)</th><th>(fs)</th><th>(sn)</th><th>(sf)</th></tr></thead><tbody><tr><td>R<sup>2</sup></td><td>0.982</td><td>0.998</td><td>0.997</td><td>0.993</td><td>0.998</td></tr></tbody></table>

The closest to 1 R<sup>2</sup> is, the best the fit of the experimental data by the theoretical data will be. The aim is thus to find the set of material parameters $(\alpha_1, \alpha_2, \alpha_3, \alpha_4, \beta_1, \beta_2, \beta_3, \beta_4)$ minimizing the ratio $S_{res}/S_{tot}$. The experimental data were extracted from [1] and the theoretical values were computed by using the stress expressions of Table 3 in combination with Eq. (29):

$$ \sigma_{nf} = \sigma_{ns} = 2k^3 \alpha_3 \exp(\beta_3 k^4) \qquad (45) $$

$$ \sigma_{fn} = 2k^3 \left[ \alpha_1 \exp(\beta_1 k^4) + 2\alpha_4 (k^2 + 2) (k^2 + 4) \exp(\beta_4 k^4 (k^2 + 4)^2) \right] \qquad (46) $$

$$ \sigma_{fs} = 2k^3 \left[ \alpha_1 \exp(\beta_1 k^4) + 2\alpha_4 (k^2 + 4) (k^2 + 8) \exp(\beta_4 k^4 (k^2 + 8)^2) \right] \qquad (47) $$

$$ \sigma_{sn} = 2k^3 \left[ \alpha_2 \exp(\beta_2 k^4) + 2\alpha_4 (k^2 + 2) (k^2 + 4) \exp(\beta_4 k^4 (k^2 + 4)^2) \right] \qquad (48) $$

$$ \sigma_{sf} = 2k^3 \left[ \alpha_2 \exp(\beta_2 k^4) + 2\alpha_4 (k^2 + 4) (k^2 + 8) \exp(\beta_4 k^4 (k^2 + 8)^2) \right] \qquad (49) $$

We have implemented equations Eqs. (45)-(49) in the Octave free software [14] and performed the material identification by using the optimization routine `fminunc` based on gradient search. The identified material parameters are shown in Table 4 and the corresponding coefficients of determination R<sup>2</sup> in Table 5. It is observed an excellent level of correlation, confirmed by Fig. 3 which shows a very good fit between the theoretical curves and the experimental data. That proves that a model using polyconvex invariants is as efficient as a model using standard invariants for predicting the mechanical response of passive ventricular myocardium. The worst prediction, with a coefficient of determination equal to 0.982, is obtained with the (nf) = (ns) modes, but it is not surprising because these modes concern the lowest level of stress and it is generally difficult for a model to match perfectly experimental low values. Finally, we have observed that the theoretical curves drawn in Fig. 3 are identical to the ones provided in Figure 6 of [2], proving the consistency of our model with the literature.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison between the model and the experimental data.</figcaption></figure>

## 7. Convexity of the SEF defined by Eq. (26)

The convexity of a SEF is an important issue regarding computational efficiency, and particularly for ensuring the uniqueness of the solution of a hyperelastic problem. In order to present sufficient conditions providing this desired property of convexity, we proceeded at the beginning in the same way as in [2]. Since the SEF defined by Eq. (26) is expressed as a sum over exponential functions, a sufficient condition for ensuring convexity is to study separately the convexity for each term. The contribution in the SEF of the first invariant $L_1$ is for example written as follows from Eq. (26):

$$W_1(C) = f_1(L_1(C)) \qquad (50)$$

With the real-valued function $f_1$ defined by:

$$f_1(x) = \frac{\alpha_1}{2\beta_1} \{ \exp[\beta_1(x-1)^2] - 1 \} \qquad (51)$$

From Eq. (8)–(10), the invariant $L_1(C)$ can be expressed linearly with respect to C thanks to the trace operator Tr:

$$L_1(C) = \langle C f_0, f_0 \rangle = \mathrm{Tr}(f_0 \otimes f_0 C) \qquad (52)$$

Deriving twice Eq. (50), one obtains from Eq. (52):

$$W_1''(C)(dC)(dC) = f_1''(L_1(C))[\mathrm{Tr}(f_0 \otimes f_0 dC)]^2 \qquad (53)$$

where the second derivative of $f_1$ is obtained from Eq. (51):

$$f_1''(L_1(C)) = \alpha_1 \{1 + 2\beta_1(L_1(C) - 1)^2\} \exp[\beta_1(L_1(C) - 1)^2] \qquad (54)$$

It is deduced from Eqs. (53) and (54) that a sufficient condition for $W_1$ to be strictly convex is:

$$\alpha_1 > 0; \quad \beta_1 > 0 \qquad (55)$$

Following the same logic with the second and third polyconvex invariants $L_2$ and $L_3$, it is obtained the additional sufficient conditions:

$$\alpha_2 > 0; \quad \alpha_3 > 0; \quad \beta_2 > 0; \quad \beta_3 > 0 \qquad (56)$$

The case of the fourth invariant is more complicated because $L_4(C)$ is expressed in a quadratic form with respect to C:

$$L_4(C) = 4[\mathrm{Tr}(f_0 \otimes s_0 C)]^2 + [\mathrm{Tr}(f_0 \otimes f_0 C) + \mathrm{Tr}(s_0 \otimes s_0 C)]^2 \qquad (57)$$

From Eq. (26), the contribution in the SEF of $L_4$ is:

$$W_4(C) = f_4(L_4(C)) \qquad (58)$$

With the real-valued function $f_4$ defined by:

$$f_4(x) = \frac{\alpha_4}{2\beta_4} \{ \exp[\beta_4(x-4)^2] - 1 \} \qquad (59)$$

Deriving twice Eqs. (57), (58) and (59), one obtains:

$$W_4''(C)(dC)(dC) = f_4''(L_4(C))[L_4'(C)(dC)]^2 + f_4'(L_4(C))L_4''(C)(dC)(dC) \qquad (60)$$

$$f_4'(L_4) = \alpha_4 (L_4 - 4) \exp[\beta_4 (L_4 - 4)^2] \qquad (61)$$

$$f_4''(L_4) = \alpha_4 (1 + 2\beta_4 (L_4 - 4)^2) \exp[\beta_4 (L_4 - 4)^2] \qquad (62)$$

$$L_4''(C)(dC)(dC) = 8[\mathrm{Tr}(f_0 \otimes s_0 dC)]^2 + 2[\mathrm{Tr}(f_0 \otimes f_0 dC) + \mathrm{Tr}(s_0 \otimes s_0 dC)]^2 \qquad (63)$$

Noting that $L_4''(C)(dC)(dC) \ge 0$, as a sum over positive quantities, and that $L_4 - 4 \ge 0$ in each of the 6 shear loadings studied in this work (last row of Table 2), it results from Eqs. (60)–(63) that a sufficient condition for $W_4$ to be convex is:

$$\alpha_4 > 0; \quad \beta_4 > 0 \qquad (64)$$

Note that the material parameters identified in Table 4 satisfy the convexity conditions given by Eqs. (55), (56) and (64).

## 8. Conclusions

In this study, a new strain energy function (SEF) was proposed for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. The SEF was built by combining exponential functions as proposed in [2] but, instead of using standard mixed invariants, we have selected the integrity basis of polyconvex invariants introduced in [10].

Using this new set of invariants allows to replace the classical mixed invariant $I_8$, which is demonstrated to be non-polyconvex, by the polyconvex invariant $L_4$ defined by Eqs. (8)–(10). With a SEF built on this basis, according to Eqs. (10) and (26), we demonstrated that our model is capable to accurately predict the experimental data obtained in [1] with 6 different shear modes applied to passive ventricular myocardium. In addition, our model is also consistent with results previously obtained in the literature [2]. Besides, we established sufficient conditions ensuring the convexity of the SEF and expressed them as inequalities involving the material parameters (Eqs. (55), (56) and (64)). Finally, our model can explicitly represent the stretch effect in the 3 specific directions of a passive ventricular myocardium: the myocyte axis direction $f_0$, the direction $s_0$ lying within the muscle layer and transverse to $f_0$, and the direction $n_0$ normal to the muscle layer. We demonstrated that 3 of the polyconvex invariants are equal to the squared stretches $I_{4f}$, $I_{4s}$ and $I_{4n}$. Consequently, three terms directly connected to the 3 possible stretch effects are embedded in the stress tensor (Eq. (34)). This is a significant difference with the model proposed in [2] which uses the 3 invariants $I_1$, $I_{4f}$ and $I_{4s}$ instead of $I_{4f}$, $I_{4s}$ and $I_{4n}$.

Because polyconvexity is often considered as a prerequisite for ensuring the existence of solutions, the next step of our research work will be to implement the new SEF introduced in this paper in a university finite element code.

### CRediT authorship contribution statement

**Renye Cai:** Methodology, Software, Formal Analysis, Investigation, Writing - review & editing.
**Frédéric Holweck:** Conceptualization, Writing - review & editing, Supervision.
**Zhi-Qiang Feng:** Conceptualization, Writing - review & editing, Supervision.
**François Peyraut:** Conceptualization, Methodology, Software, Formal analysis, Writing - original draft, Writing - review & editing, Visualization, Supervision, Project administration.
Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

References

[1] S. Dokos, B.H. Smaill, A.A. Young, I.J. LeGrice, Shear properties of pas-
sive ventricular myocardium, Am. J. Physiol. Heart Circ. Physiol. (2002)
283H2650–H2659.

[2] G.A. Holzapfel, R.W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Phil. Trans. R. Soc. A 367 (2009) 3445–3475.

[3] T.S.E. Eriksson, A.J. Prassl, G. Plank, G.A. Holzapfel, Modeling the dispersion in electromechanically coupled myocardium, Int. J. Numer. Methods Biomed. Eng. 29 (2013) 1267–1284.

[4] A.V. Melnik, X. Luo, R.W. Ogden, A generalised structure tensor model for the mixed invariant $I_s$, Int. J. Non-Linear Mech. 107 (2018) 137–148.

[5] E. McEvoy, G.A. Holzapfel, P. McGarry, Compressibility and anisotropy of the ventricular myocardium: Experimental analysis and microstructural modeling, J. Biomech. Eng. 140 (2018) 081004-1.

[6] A. Palit, S.K. Bhudia, T.N. Arvanitis, G.A. Turley, M.A. Williams, In vivo estimation of passive biomechanical properties of human myocardium, Med. Biol. Eng. Comput. 56 (2018) 1615–1631.

[7] D.S. Li, R. Avazmohammadi, S.S. Merchant, T. Kawamura, E.W. Hsu, J.H. Gorman III, R.C. Gorman, M.S. Sacks, Insights into the passive mechanical behavior of left ventricular myocardium using a robust constitutive model based on full 3D kinematics, J. Mech. Behav. Biomed. Mater. 103 (2020) 103508.

[8] D.J. O'Shea, M.M. Attard, D.C. Kellermann, Hyperelastic constitutive modelling for transversely isotropic composites and orthotropic biological tissues, Int. J. Solids Struct. 169 (2019) 1–20.

[9] A.T. Ta, N. Labed, F. Holweck, A. Thionnet, F. Peyraut, A new invariant-based method for building biomechanical behavior laws–application to an anisotropic hyperelastic material with two fiber families, Int. J. Solids Struct. 50 (14) (2013) 2251–2258.

[10] R. Cai, F. Holweck, Z.Q. Feng, F. Peyraut, A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials, Int. J. Solids Struct. 115–116 (2017) 126–139.

[11] A. Thionnet, C. Martin, A new constructive method using the theory of invariants to obtain material behavior laws, Int. J. Solids Struct. 43 (2) (2006) 325–345.

[12] J.M. Ball, Convexity conditions and existence theorems in nonlinear elasticity, Arch. Ration. Mech. Anal. 63 (4) (1976) 337–403.

[13] J. Schröder, P. Neff, Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions, Int. J. Solids Struct. 40 (2) (2003) 401–445.

[14] J.W. Eaton, D. Bateman, S. Hauberg, R. Wehbring, GNU Octave version 4.2.0 manual: a high-level interactive language for numerical computations, 2016.