A modification of Holzapfel-Ogden hyperelastic model of myocardium better describing its passive mechanical behavior

Contents lists available at ScienceDirect European Journal of Mechanics / A Solids journal homepage: www.elsevier.com/locate/ejmsol Full length article A modification of Holzapfel–Ogden hyperelastic model of myocardium better describing its passive mechanical behavior Jiří Vaverka ∗, Jiří Burša Institute of Solid Mechanics, Mechatronics and Biomechanics, Faculty of Mechanical Engineering, Brno University of Technology, Brno, Czech Republic ARTICLE INFO Keywords: Cardiac mechanics Myocardium Hyperelasticity Constitutive model Orthotropy ABSTRACT The passive mechanical behavior of the myocardium is usually mathematically described within the framework of hyperelasticity. One of the most popular models of this kind is that proposed by Holzapfel and Ogden in 2009. It is an orthotropic model formulated in terms of a reasonably selected set of scalar invariants representing different components of the myocardium. Several modifications of the model have emerged over the years. In this paper, we present another one which is characterized by an innovative approach to the modeling of myocardial ‘‘sheets’’, i.e. lamellar collagenous structures that endow the myocardium with orthotropic mechanical properties. We describe their contribution by means of a less common scalar invariant which expresses the change of area of an oriented planar element (representing the plane of a sheet). To compare our formulation with the original model, we matched both of them to the biaxial tension and simple shear experimental data from the literature using a nonlinear least-squares optimization algorithm. The objective function for each model included both biaxial and simple shear data in order to obtain a single set of parameters for both deformation modes. The results show that our modified model can accurately describe both types of tests. The total residual is lowered by approximately 80% by our modification and 𝑅2increases from 0.877 to 0.978 which demonstrates the significance of our modification on the quality of the fit. 1. Introduction Accurate mathematical description of the mechanical behavior of passive (relaxed) myocardium is essential in order to study mechanical phenomena in the heart by computational modeling. Constitutive equations providing such description are particularly important for the development of computer models based on the finite element method. Computational studies published over the last few decades proved that such models can substantially increase our knowledge and understanding of the tissue properties and physiological processes in a healthy heart (e.g., Xi et al.,2019;Kerckhoffs et al.,2003;McEvoy et al.,2018), help us to explain mechanisms underlying heart diseases and investigate their consequences (e.g., Kovacheva et al.,2021;Mojumder et al., 2023;Costabal et al.,2019), as well as allow us to predict the outcomes of medical treatments and interventions (e.g., Xu et al.,2020;Li et al., 2020a;Dabiri et al.,2019). The vast majority of these works treat myocardium as a hyperelastic material despite the reported viscoelastic features (rate-dependence, hysteresis) exhibited by myocardium during mechanical tests (see, e.g., Sommer et al.,2015;Dokos et al.,2002). However, a recent computational study by Tikenoğulları et al. (2022) has justified this simplification by showing that viscous relaxation of ∗Correspondence to: Institute of Solid Mechanics, Mechatronics and Biomechanics, Brno University of Technology, Technická 2896/2, 616 69, Brno, Czech Republic. E-mail address: [email protected] (J. Vaverka). the myocardium has a negligible effect on the overall behavior of the whole heart during physiologically relevant time scales of the cardiac cycle. Thus the approximation by a hyperelastic model is acceptable for most practical applications, provided the chosen model can accurately capture the orthotropic and highly nonlinear (elastic) response of myocardium revealed by experiments (Sommer et al.,2015;Dokos et al.,2002). Considering further that nonlinear viscoelastic models of the myocardium are highly complex and include additional material parameters that are difficult to estimate (see, e.g., Gültekin et al.,2016; Zhang et al.,2023), it can be expected that hyperelastic models will remain dominant in the field of computational cardiac biomechanics within the near future; thus their further development is still important. Obviously, a hyperelastic model is more likely to provide an accurate mechanical description of myocardium if the strain–energy function (SEF) defining the model reflects the internal microstructure of the tissue. This can be efficiently achieved by formulating the function as a sum of separate terms, each of which depends only on a single scalar invariant with a clear physical interpretation. Each term of the model then represents certain constituent(s) of the tissue and the character of https://doi.org/10.1016/j.euromechsol.2025.105586 Received 18 August 2024; Received in revised form 30 December 2024; Accepted 19 January 2025 European Journal of Mechanics / A Solids 111 (2025) 105586 Available online 27 January 2025 0997-7538/© 2025 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). J. Vaverka and J. Burša its dependence on the corresponding invariant (e.g., quadratic or exponential) is determined on the basis of experimental stress–strain data. Such a rational, structurally and experimentally motivated approach was adopted also by Holzapfel and Ogden (2009) who proposed an orthotropic hyperelastic model for myocardium (hereinafter abbreviated as the “HO” model) which has since then become one of the most popular models in the field. However, several modified versions of the model have emerged over the years and some of them are preferred nowadays by many researchers over the original formulation. One of the first (and most notable) modifications was proposed by Göktepe et al. (2011) who introduced the traditional isochoric–volumetric decoupling (see Holzapfel,2000) into the (originally incompressible) HO model. This was effected by replacing all invariants in the SEF with their isochoric versions and extending the function with a volumetric term. The resulting nearly-incompressible formulation was then used by many other authors (e.g., Eriksson et al.,2013a;Palit et al.,2018); besides, it was also implemented in commercial finite element software Abaqus®. In some other works the isochoric invariants were introduced only to the isotropic part of the SEF while the anisotropic part was left unchanged (e.g., Gerbi et al.,2019;McEvoy et al.,2018); this was motivated by the works of Nolan et al. (2014) and Vergori et al. (2013) who have shown that isochoric anisotropic invariants produce unphysical responses under certain deformation modes. An important improvement was made by Eriksson et al. (2013b) and Melnik et al. (2018) who introduced fiber dispersion into anisotropic invariants. Other significant modifications include a substitution of the exponential isotropic term by a polynomial one (McEvoy et al.,2018), simplification of the model by omitting two of its three anisotropic terms (Chapelle and Le Gall,2023) or, in contrast, its extension by two more anisotropic terms (with four additional material parameters) (Li et al.,2020b). The most general HO model with 6 anisotropic terms (instead of 3) was considered by Guan et al. (2019) who studied its ability to describe three different experimental data sets from literature. By analyzing the contribution of each term to the overall goodness of fit, they proposed three reduced versions of the model (with 3–4 anisotropic terms) sufficient for an accurate description of the three data sets. The same general model was investigated also by Martonová et al. (2024). The above examples show that many different versions of the HO model can already be found in the literature, each of them designed to fit the specific needs of its authors and, therefore, more suitable for particular applications than the others. In this paper, we present a new modification which is distinctive from any other formulation presented in the literature by the manner in which it deals with the fact that the ventricular muscle fibers are aggregated into the so-called “sheets” (e.g., Stephenson et al.,2018) that provide the myocardium with orthotropic rather than transversely isotropic mechanical properties. Unlike Holzapfel and Ogden (2009), we describe the mechanical contribution of sheets by an exponential strain–energy term proposed by Balzani et al. (2006) which is formulated in terms of a (not very usual) scalar invariant 𝐾1which is a measure of change in relation to a deformation of the area of a preferred planar element. We consider this invariant to be a suitable measure of the deformation of sheets, similar in its nature to the analogous invariants that express stretch of fibers or change of volume and that are regularly used in constitutive equations. Since we introduce no other change to the model, our modification retains the mathematical structure of the HO model, its exponential form, as well as the number of material parameters. However, we will show that it has considerably improved fitting capabilities compared to the original formulation. We demonstrate this by fitting the model to data from five different biaxial tension tests and six different simple shear tests (all taken from Sommer et al.,2015). The comparison with the original HO model reveals that our modification enables a description of all the data with a single set of material parameters, which the original formulation cannot satisfactorily achieve. The paper is organized as follows. Section 2briefly summarizes the most important anatomical features of the ventricular myocardium and describes how Holzapfel and Ogden (2009) incorporated them into their successful model. Section 3presents our approach to the constitutive modeling of the myocardium, characterized by the employment of invariant 𝐾1whose physical meaning we explain in more detail. We introduce a modified SEF for the myocardium and derive the corresponding relation for the Cauchy stress. Following Holzapfel and Ogden (2009), we assume fully incompressible behavior. In Section 4we first describe the biaxial and simple shear data published by Sommer et al. (2015), which are suitable for the evaluation of material parameters of orthotropic models, and then we discuss the controversy surrounding the correct interpretation of the cross-fiber response obtained from biaxial tests. In Section 5we fit our model and the HO model to the data from Sommer et al. (2015) in order to demonstrate the benefits of our approach. Finally, in Section 6, we discuss some problematic points, make further suggestions and formulate conclusions. 2. Anatomy of ventricular myocardium and the model by Holzapfel & Ogden Myocardium is composed mostly of cardiomyocytes, i.e. elongated contractile cells embedded in an extracellular matrix of connective tissue (Smerup et al.,2009). Cardiomyocytes are linked end-to-end into longitudinal chains, or fibers, which are interconnected by sidebranches (Smerup et al.,2009). The fibers join and split along their paths and have neither a discernible beginning nor ending (Lunkenheimer and Niederer,2012;Stephenson et al.,2018); nevertheless, in a typical small volume of the network, it is possible to distinguish the predominant fiber direction, which can be represented by a unit vector 𝐟(see Fig. 1). Adjacent fibers also have a higher degree of organization since they are laterally aggregated (through endomysial connective tissue) into larger blocks with lamellar appearance (Lunkenheimer and Niederer,2012;LeGrice et al.,1995;Stephenson et al.,2018). These higher-order structures, which are frequently called “sheets”, provide the myocardium with orthotropic mechanical properties (Sommer et al.,2015;Dokos et al.,2002) rather than transversely isotropic. The structural integrity of layers is provided mainly by endomysial collagen fibers which act as lateral connections between muscle fibers and thus enable transmission of forces in axial as well as in transverse directions (Lunkenheimer and Niederer,2012;Weber,1989). Adjacent sheets are separated by elongated perimysial spaces (also called perimysial clefts or cleavage planes) filled with gelatinous lubricating fluid which is supposed to facilitate shearing of sheets relative to each other (Lunkenheimer and Niederer,2012;Weber,1989). However, in order to prevent an excessive slippage of layers or even a rupture of the tissue, the cleavage planes are spanned by a sparse network of perimysial collagen fibers which connect adjacent layers and hence strengthen the structure (Weber et al.,1994;Weber,1989, see also micrographs in Lunkenheimer and Niederer,2012 and LeGrice et al., 1995). The orientation of a sheet in the myocardium can be mathematically described by choosing one of the two unit vectors parallel to the plane of the sheet and perpendicular to the 𝐟direction. The chosen vector is then referred to as the sheet vector, 𝐬. The directions of material orthotropy are completed with the sheet-normal vector, 𝐧, which is the vector for which the triple (𝐟,𝐬,𝐧)forms a right-handed orthonormal basis. Respecting the structural organization of the myocardium described above, Holzapfel and Ogden (2009) described its mechanical behavior in terms of invariants 𝐼1∶=  𝐼1(𝐂) ∶= tr𝐂,(1) 𝐼4f∶=  𝐼4f(𝐂) ∶= 𝐟⋅(𝐂𝐟),(2) 𝐼4s∶=  𝐼4s(𝐂) ∶= 𝐬⋅(𝐂𝐬),(3) 𝐼8∶=  𝐼8(𝐂) ∶= 𝐟⋅(𝐂𝐬),(4) European Journal of Mechanics / A Solids 111 (2025) 105586 2 J. Vaverka and J. Burša Fig. 1. Schematic representation of the lamellar structure of ventricular myocardium with corresponding structural vectors 𝐟(fiber), 𝐬(sheet), and 𝐧(normal to sheet). Created on the basis of scanning electron microscopic images published by Lunkenheimer and Niederer (2012) and LeGrice et al. (1995). which depend on the deformation gradient 𝐅through the right Cauchy– Green tensor 𝐂∶= 𝐅⊤𝐅. (Note that in (1)−(4), and in some other equations below, we use the circumflex “ˆ” to distinguish a function from its value.) As the authors explained, 𝐼1was employed to describe the deformation of all non-collagenous and non-muscular constituents (like elastin and fluids) which can reasonably be considered isotropic. The more pronounced stiffness exhibited by the myocardium in the 𝐟direction was assumed to be due to the muscle fibers as well as to the collagen network (both endomysial and perimysial). The combined effect of these two constituents was described by 𝐼4f, which expresses the square of stretch in the 𝐟direction since 𝐼4f=𝐟⋅(𝐂𝐟) =𝐟⋅(𝐅⊤𝐅𝐟) = (𝐅𝐟)⋅(𝐅𝐟) =‖𝐅𝐟‖2.(5) Additional stiffness in the 𝐬direction, associated with the collagen fibers connecting the muscle fibers, was modeled through 𝐼4swhose physical meaning is, of course, analogous to 𝐼4f. Finally, 𝐼8was included in order to stiffen the response to the simple shear modes (fs) and (sf), and thereby to distinguish them from the modes (fn) and (sn), see Fig. 2. The decision of Holzapfel and Ogden (2009) to increase stiffness in modes (fs) and (sf) was motivated by the experimental data of Dokos et al. (2002) who tested cubic samples from pig myocardium in all the six possible simple shear modes; they reported considerably stiffer behavior of the myocardium in the (fs) mode compared to the (fn) mode, and in the (sf) mode compared to the (sn) mode. Although the results from the simple shear tests of human myocardium, published later by Sommer et al. (2015), did not confirm such significant differences, no complete data set comparable to that of Dokos et al. (2002) was available at the time of publication of the HO model, therefore the inclusion of the 𝐼8invariant seemed necessary. The sensitivity of 𝐼8to the shear modes (fs) and (sf) becomes evident if we rewrite (4) as 𝐼8=𝐟⋅(𝐂𝐬) =𝐟⋅(𝐅⊤𝐅𝐬) = (𝐅𝐟)⋅(𝐅𝐬) = cos (𝜃)‖𝐅𝐟‖‖𝐅𝐬‖,(6) where 𝜃is the angle between 𝐅𝐟 and 𝐅𝐬. Thus 𝐼8= 0in all simple shear modes except for (fs) and (sf). In order to reproduce exponential trends shown by experimental data, Holzapfel and Ogden (2009) embedded all four invariants into exponential functions based on that proposed by Demiray (1972). The resulting SEF 𝜓HO of the argument HO ∶=  HO(𝐂) ∶= (𝐼1, 𝐼4f, 𝐼4s, 𝐼8)(7) is given by 𝜓HO(HO) ∶= 𝑎1 2𝑏1 (exp(𝑏1(𝐼1− 3)) − 1) +𝑎f 2𝑏f (exp(𝑏f(𝐼4f− 1)2) − 1) +𝑎s 2𝑏s (exp(𝑏s(𝐼4s− 1)2) − 1) +𝑎fs 2𝑏fs (exp(𝑏fs𝐼2 8) − 1), (8) where 𝑎1,𝑎f,𝑎sand 𝑎fs are stress-like material parameters, while 𝑏1, 𝑏f,𝑏sand 𝑏fs are dimensionless. To ensure the existence of minimizers in boundary value problems and to obtain physically meaningful responses, the 𝑎-parameters must be non-negative and the 𝑏-parameters must be positive (see the discussion in Section 6 of Holzapfel and Ogden,2009). Also, each of the middle two terms in (8) should be set to zero when its corresponding invariant (𝐼4for 𝐼4s) is less than 1 which ensures convexity of both terms in the full range of deformations (Balzani et al.,2006). This deactivation is also physically reasonable because muscle and collagen fibers cannot be expected to support significant compression. Identification of the parameters of the HO model on the basis of experimental data has been performed by several authors. However, most of them either used only simple shear data (e.g., Göktepe et al., 2011;Wang et al.,2013;McEvoy et al.,2018), or they used both simple shear data and biaxial data but they provided different parameters for each test type (e.g., Holzapfel and Ogden,2009;Gültekin et al.,2016). This can be attributed to the fact that a fusion of biaxial and simple shear data in one fitting process does not give satisfactory results with the HO model, as we will show in Section 5. This deficiency of the model was recognized earlier by Gültekin et al. (2016) (although in the context of viscoelasticity) and it is also apparent from the fitting results of Guan et al. (2019) and Martonová et al. (2024). Our approach presented in this paper aims to overcome this deficiency of the HO model without increasing the number of parameters. 3. The proposed model It is reasonable to expect that a typical sheet (composed of endomysial tissue and muscle fibers) resists forces acting in any direction within its plane. This assumption is consistent with Holzapfel and Ogden (2009) who assumed that collagen fibers increase the stiffness in both 𝐟and 𝐬directions. The mechanical properties of endomysial connective tissue might, of course, be anisotropic (because of the potentially nonuniform directional distribution of collagen fibers) but there seems to be no relevant information on this issue; thus we find it appropriate to idealize a typical sheet as an isotropic plane reinforced by a family of parallel muscle fibers. This planar structure can then be imagined as being embedded in a (presumably) isotropic three-dimensional matrix representing the rest of the extracellular components (including the perimysial tissue). This modeling approach leaves us with three idealized structural components of the myocardium whose mechanical behavior should now be described in terms of suitable invariants inserted into separate, additively contributing strain–energy terms. We see no good reason to change anything on the description of extracellular matrix and muscle fibers proposed by Holzapfel and Ogden (2009); thus we retain the first two terms in Eq. (8) but we associate the 𝐼4fterm solely with muscle fibers and not with collagen fibers. However, we take a different approach to the description of deformation and mechanical response of collagenous endomysium, which is based on a not very common invariant which Schröder and Neff (2003) and Balzani et al. (2006) denote as 𝐾1. Using the definition of the cofactor of any invertible tensor 𝐋, cof(𝐋) ∶= det(𝐋)𝐋−⊤,(9) the invariant is defined by 𝐾1∶=  𝐾1(𝐂) ∶= 𝐧⋅(cof(𝐂)𝐧).(10) Vector 𝐧in (10) could be an arbitrary vector determining a preferred material direction, but here we will regard it as representing the sheetnormal direction illustrated in Fig. 1. This choice will give 𝐾1a physical meaning suitable for our purpose, which can be disclosed using the following result (cf. Eq. (5)): 𝐾1=𝐧⋅(cof(𝐂)𝐧) =𝐧⋅(cof(𝐅)⊤cof(𝐅)𝐧) = (cof(𝐅)𝐧)⋅(cof(𝐅)𝐧) =‖cof(𝐅)𝐧‖2.(11) European Journal of Mechanics / A Solids 111 (2025) 105586 3 J. Vaverka and J. Burša Fig. 2. Six possible modes of simple shear deformation of a cubic myocardial specimen relative to the principal material directions 𝐟,𝐬and 𝐧. We specify a simple shear mode by a parenthetical symbol in which the first letter refers to the normal to the face of the cubic specimen that is shifted by the simple shear and the second denotes the direction of shift. Dashed lines represent cube edges in the reference configuration. The figure was drawn on the basis of similar illustrations used elsewhere; e.g., in Holzapfel and Ogden (2009). Fig. 3. Illustration of the action of 𝐅and cof(𝐅). Deformation of a neighborhood of a material point 𝐗onto a neighborhood of the corresponding point 𝐱transforms material vectors 𝐟,𝐬and 𝐧respectively to 𝐅𝐟,𝐅𝐬 and 𝐅𝐧. However, cof(𝐅)maps 𝐧(regarded as area vector) to cof(𝐅)𝐧which is perpendicular to the parallelogram with area  𝑃 (initially 𝑃) and its magnitude equals  𝑃. Since 𝐧=𝐟×𝐬, we can regard 𝐧as the area vector of the parallelogram defined by 𝐟and 𝐬. Then, using the identity cof(𝐅)(𝐟×𝐬) = (𝐅𝐟) × (𝐅𝐬)(12) (cf., e.g., Gurtin et al.,2010, p. 23), we can see that cof(𝐅)𝐧is the area vector of the parallelogram defined by 𝐅𝐟 and 𝐅𝐬 (i.e., the vector orthogonal to the parallelogram, directed according to the right-hand screw rule, and with magnitude equal to the area of the parallelogram; cf. Fig. 3). Thus, by (11),𝐾1is the square of the area of the parallelogram, or, since ‖𝐧‖= 1, the square of the relative change of the area of a planar element initially perpendicular to 𝐧. We use 𝐾1in a function having the same exponential form as the terms with 𝐼4fand 𝐼4sin the HO model (8). The usage of this term (the third one in Eq. (14)) in SEFs for soft biological tissues was suggested by Balzani et al. (2006) but we are not aware of any previous attempt to employ it for the description of the myocardium. Since the term can capture the collagen-induced stiffening in the 𝐬direction, there is no further need for the 𝐼4sinvariant and thus we omit the corresponding term of the HO model. However, in order to make the model adaptable to qualitatively different simple shear data available in the literature (cf. Section 2), we maintain the last term of the HO model containing the 𝐼8invariant. The resulting modified model is then defined in terms of the family of four invariants ∶=  (𝐂) ∶= (𝐼1, 𝐼4f, 𝐾1, 𝐼8)(13) by a SEF 𝜓() ∶= 𝑎1 2𝑏1 (exp(𝑏1(𝐼1− 3)) − 1) +𝑎f 2𝑏f (exp(𝑏f(𝐼4f− 1)2) − 1) +𝑎nn 2𝑏nn (exp(𝑏nn(𝐾1− 1)2) − 1) +𝑎fs 2𝑏fs (exp(𝑏fs𝐼2 8) − 1). (14) Material parameters 𝑎1,𝑎f,𝑎nn and 𝑎fs are stress-like, while 𝑏1,𝑏f, 𝑏nn and 𝑏fs are unitless. In order to satisfy the requirements related to convexity, the 𝑎-parameters must be non-negative, the 𝑏-parameters must be positive, and each of the middle two terms in (14) should be set to zero when its corresponding invariant (𝐼4for 𝐾1) is less than 1 (Balzani et al.,2006). For clarity, we will denote by  𝛹the SEF 𝜓 expressed as a function of 𝐂instead of , i.e.:  𝛹∶= 𝜓◦ .(15) Then, with the assumption of incompressibility, the second Piola– Kirchhoff stress 𝐒corresponding to 𝜓 can be expressed as 𝐒∶=  𝐒(𝐂) ∶= 2∇𝐂 𝛹−𝑝𝐽𝐂−1 (16) and the Cauchy stress 𝝈is given by 𝝈∶= 𝐽−1𝐅𝐒𝐅⊤= 2𝐽−1𝐅(∇𝐂 𝛹)𝐅⊤−𝑝𝟏.(17) The presence of the volume change 𝐽∶= det(𝐅)in (16) and (17) may seem unnecessary since 𝐽= 1for incompressible materials, but maintaining 𝐽is important for the correct evaluation of elasticity tensors (cf. Section 6.5.1 of Bonet and Wood,2008). For this reason, we keep 𝐽(or equivalently det(𝐂)) in several equations in this paper. The scalar 𝑝in (16) and (17) is an indeterminate pressure (meaning that it cannot be determined from a constitutive equation). Note that in this particular case, it does not coincide with the pressure part of 𝝈, i.e., 𝑝≠−1 3tr𝝈(cf. Section 6.5.1 of Bonet and Wood,2008). In the following, we will use the abbreviations ∶= (𝐂−1𝐧)⊗(𝐂−1𝐧),(18) European Journal of Mechanics / A Solids 111 (2025) 105586 4 J. Vaverka and J. Burša ∶= (cof(𝐅)𝐧)⊗(cof(𝐅)𝐧).(19) Using the chain rule, we can express ∇𝐂 𝛹= ∇𝜓∇𝐂  =𝜓,𝐼1()∇𝐂 𝐼1+𝜓,𝐼4f()∇𝐂 𝐼4f+𝜓,𝐾1()∇𝐂 𝐾1+𝜓,𝐼8()∇𝐂 𝐼8.(20) Gradients in the last line of (20) satisfy the identities ∇𝐂 𝐼1=𝟏,(21) ∇𝐂 𝐼4f=𝐟⊗𝐟,(22) ∇𝐂 𝐾1=𝐾1𝐂−1 −det(𝐂),(23) ∇𝐂 𝐼8= (1∕2)(𝐟⊗𝐬+𝐬⊗𝐟),(24) and partial derivatives are given by 𝜓,𝐼1() = (𝑎1∕2) exp(𝑏1(𝐼1− 3)),(25) 𝜓,𝐼4f() =𝑎f(𝐼4f− 1) exp(𝑏f(𝐼4f− 1)2),(26) 𝜓,𝐾1() =𝑎nn(𝐾1− 1) exp(𝑏nn(𝐾1− 1)2),(27) 𝜓,𝐼8() =𝑎fs𝐼8exp(𝑏fs𝐼2 8).(28) Eqs. (21)–(24) can now be substituted into (20) and the resulting expression can be used in (17) to obtain the Cauchy stress relation in an extended form 𝝈= 2𝜓,𝐼1𝐁+ 2𝜓,𝐼4f 𝐟⊗ 𝐟+ 2𝜓,𝐾1(𝐾1𝟏−) +𝜓,𝐼8( 𝐟⊗ 𝐬+ 𝐬⊗ 𝐟) −𝑝𝟏,(29) where  𝐟∶= 𝐅𝐟 and  𝐬∶= 𝐅𝐬. Note that we omitted the argument from the partial derivatives in (29) in order to avoid clutter. The general expression (29) can be used, e.g., to derive explicit stress–strain relations for deformation modes to which specimens are commonly subjected during mechanical tests. Such relations can then be matched with corresponding experimental data in order to evaluate the material parameters of the model, which is what we do in Section 5. 4. Experimental data for fitting, interpretation of the cross-fiber results of biaxial tests In order to compare our modification 𝜓, Eq. (14), with the original SEF 𝜓HO, Eq. (8), we will fit their corresponding stress relations to the data extracted from Figs. 9 and 13 in Sommer et al. (2015) which contain averaged results of five different biaxial tests (Fig. 9) and six different simple shear tests (Fig. 13) of human myocardium. The biaxial tests were performed on thin squared specimens with dimensions of 25 × 25 mm and a thickness of approximately 2.3mm. They were cut parallel to the ventricular wall so that one loading direction of each specimen was aligned with the predominant 𝐟direction and the other one was therefore approximately perpendicular to the fibers. The former direction was called the mean-fiber direction (MFD) and the latter the cross-fiber direction (CFD). Fig. 9 of Sommer et al. (2015) contains biaxial responses for five different strain ratios between the MFD and the CFD, namely: 0.5∶1,0.75∶1,1∶1 (equibiaxial test), 1∶0.75 and 1∶0.5. The maximum applied strain was 10% in all cases. The simple shear tests captured in Fig. 13 of Sommer et al. (2015) were performed on small cubic specimens excised from regions adjacent to the biaxially tested specimens in order to ensure uniform material properties. The edges of each specimen were approximately 4mm long and they were oriented along the material directions 𝐟,𝐬and 𝐧. The specimens were tested in all six possible simple shear modes associated with these directions (cf. Fig. 2). The maximum applied amount of shear was 0.5 in all modes. Before using the above-described experimental data to obtain the material parameters of the investigated models, we must first establish how to treat the responses in the CFD obtained from biaxial tests. Specifically, we must decide whether it is more appropriate to consider these results as representing the mechanical behavior of the myocardium in the 𝐬direction or in the 𝐧direction. As we discuss below, this is not a straightforward task because the orientation of sheets inside the wall is inhomogeneous and highly variable and therefore it is difficult to make a conclusive assessment of the relation between the CFD and the average 𝐬or 𝐧direction in a typical biaxial specimen. As a consequence, some authors in the past associated the CFD with the 𝐬direction (e.g., Holzapfel and Ogden,2009;McEvoy et al.,2018) while others with the 𝐧direction (e.g., Gültekin et al.,2016;Guan et al., 2019). We will now provide a short summary of relevant information from the literature about the sheet orientation that has influenced our decision on this issue. The orientation of sheets in the left ventricle is heterogeneous, as is evident, e.g., from the results of Agger et al. (2017) who analyzed in detail the myocardial architecture in ventricles by means of the diffusion tensor magnetic resonance imaging technique. Similar findings were reported also by Dokos et al. (2002) who studied transmural segments of the left ventricular wall before they cut them into cubic simple-shear specimens. They wrote: “In all hearts tested, midwall muscle layers were inclined at ∼45◦to the radial direction. Although it was possible to identify a predominant midwall orientation of myocardial laminae on the transmural cut surfaces, there were commonly regions where layers were oriented at ∼90◦to this principal direction.” Besides confirming the heterogeneous orientation of sheets, this description asserts that they are either approximately parallel to the local epicardial tangential plane, or inclined from it by approximately 45◦or less, but they are unlikely to extend in a radial fashion. Contrary to this, Sommer et al. (2015) in their Fig. 3marked the sheet directions in myocardial samples by arrows directed almost radially, but unfortunately, they did not formulate any general summary comparable to that of Dokos et al. (2002). Also, since they tested specimens from the outer, middle and inner portion of the wall, it is likely that the predominant orientation considerably varied between specimens (because of the heterogeneity). Another important finding about the orientation of sheets is that the angle between 𝐬and the local wall tangent plane significantly changes during the cardiac cycle. In particular, Ferreira et al. (2014) conducted in-vivo diffusion tensor magnetic resonance imaging of healthy human hearts and observed the mean angle of 24.0◦at end-diastole and of 56.4◦ at end-systole. A similar study was performed by Nielles-Vallespin et al. (2017) who reported the mean values of 18◦and 65◦at end-diastole and end-systole, respectively. Finally, the above-mentioned study by Agger et al. (2017) used the same technique to measure the angle between 𝐧 (rather than 𝐬) and the tangent plane in excised pig hearts fixed at the end-diastolic state. Their angle histograms show that the angle values were most frequently between 60◦and 90◦(cf. their Fig. 3) which is consistent with the low end-diastolic values (cited above) of the angle between 𝐬and the tangential plane. All these results suggest that the answer to the question of whether a biaxial test specimen is more likely to be spanned by local 𝐟and 𝐬 vectors or by 𝐟and 𝐧vectors depends on the actual configuration of the myocardial sample at the moment when the specimen is being cut from it. In this regard, it should be recalled that experimenters always do what they can to keep their tissue samples in a non-contracted state (see, e.g., the description of the preparation process in Sommer et al., 2015). Consequently, the samples are always in a relaxed state. If we add to this the fact that the myocardium is fully contracted at the end of systole, it follows that the tissue samples for experiments should be much closer to the end-diastolic state of the myocardium, which is characterized by the tangential orientation of sheets rather than radial. The same result can be obtained by comparing the left ventricle free wall thicknesses of the 28 hearts used by Sommer et al. (2015) (cf. their Table 1) with normal average thickness values at end-diastole and endsystole. Since there is a significant correlation between wall thickness changes and sheet angle changes during systole (Nielles-Vallespin et al., 2017), it is possible to infer the predominant orientation of sheets in the samples from their thicknesses. According to Peshock et al. (1989) and Dawson et al. (2011), the mean end-diastolic free wall thickness in a normal healthy heart is approximately 10 mm and the percentual systolic wall thickening reaches 50 − 60%. These reference values can be European Journal of Mechanics / A Solids 111 (2025) 105586 5 J. Vaverka and J. Burša compared with the mean wall thickness of the hearts used by Sommer et al. (2015), which is 12.5mm. If we add the fact that the majority of the patients whose hearts were tested had some record of heart-related disease(s) (including in particular 5 hearts with diagnosed hypertrophy and wall thickness as high as 19 mm), it can be concluded that the wall samples were in a state which was much closer to the end-diastolic configuration of the ventricle than to the end-systolic one. The above facts and considerations lead us to believe that the stress responses obtained from biaxial tests represent the behavior of the 𝐟 𝐬 plane of myocardium rather than the 𝐟 𝐧plane, or at least that the responses in the CFD are indeed significantly influenced by the stiffness of collagenous sheets. Therefore, we would find it inappropriate to neglect the mechanical contribution of sheets during biaxial extension by assuming that they are perpendicular to the CFD (and thus also perpendicular to the heart wall which, as the above-reported angle measurements clearly show, virtually never occurs in a healthy heart). For these reasons, we decided to match the biaxial results in the CFD with the myocardial 𝐬direction and in the next section we will fit both 𝜓 and 𝜓HO using this assumption. 5. Fitting to experimental data We utilized both biaxial and simple shear data in a single optimization problem aimed at minimizing the objective function defined as the sum of weighted squares of differences between the measured and model-predicted stress values. The stress responses for our SEF 𝜓 were calculated from the general stress relation (29). The biaxial stress components corresponding to the applied stretches 𝜆fand 𝜆sare given by 𝜎f= 2((𝜆2 f−𝜆−2 f𝜆−2 s)𝜓,𝐼1+𝜆2 f𝜓,𝐼4f+𝜆2 f𝜆2 s𝜓,𝐾1),(30) 𝜎s= 2((𝜆2 s−𝜆−2 f𝜆−2 s)𝜓,𝐼1+𝜆2 f𝜆2 s𝜓,𝐾1),(31) while the shear responses corresponding to the amount of shear 𝛾are given by (fs) ∶𝜏= 2𝛾(𝜓,𝐼1+𝜓,𝐼4f) +𝜓,𝐼8,(32) (fn) ∶𝜏= 2𝛾(𝜓,𝐼1+𝜓,𝐼4f+𝜓,𝐾1),(33) (sf) ∶𝜏= 2𝛾 𝜓,𝐼1+𝜓,𝐼8,(34) (sn) ∶𝜏= 2𝛾(𝜓,𝐼1+𝜓,𝐾1),(35) (nf) ∶𝜏= 2𝛾 𝜓,𝐼1,(36) (ns) ∶𝜏= 2𝛾 𝜓,𝐼1.(37) Analogous expressions for the HO model can be found in the original paper (Holzapfel and Ogden,2009). From each of the six available simple shear responses, we extracted 𝑛p∶= 21 data points, evenly spaced in the tested range [0,0.5] of the amount of shear. Thus, in general, for the 𝑗th simple shear mode, where 𝑗∈ {1,…,6}, we obtained a family ((𝛾(𝑗) 𝑖, 𝜏(𝑗) 𝑖) ∣𝑖∈ {1,…, 𝑛p}) of 𝑛ppairs (𝛾(𝑗) 𝑖, 𝜏(𝑗) 𝑖), each consisting of the amount-of-shear value, 𝛾(𝑗) 𝑖, and the corresponding shear stress value, 𝜏(𝑗) 𝑖. (To avoid confusion, we note that we use 𝑗∈ {1,…,6} merely as a distinguishing label, the enumeration of the shear modes may be arbitrary. Notice also that we use a superposed bar to distinguish experimental stresses from the predicted ones.) Similarly, for each of the five available biaxial test records, we extracted 𝑛pdata points from the response in the MFD and the same number of points from the response in the CFD (which we identify with the 𝐬direction). In this manner, for the 𝑘th biaxial test, with 𝑘∈ {1,…,5}, we obtained a family (((𝜀f)(𝑘) 𝑖,(𝜎f)(𝑘) 𝑖) ∣𝑖∈ {1,…, 𝑛p}) of 𝑛ppairs ((𝜀f)(𝑘) 𝑖,(𝜎f)(𝑘) 𝑖), each consisting of the value of engineering strain in the 𝐟direction, (𝜀f)(𝑘) 𝑖, and the corresponding value of Cauchy stress, (𝜎f)(𝑘) 𝑖(the order of tests is again immaterial). Analogous data sets of the form (((𝜀s)(𝑘) 𝑖,(𝜎s)(𝑘) 𝑖) ∣𝑖∈ {1,…, 𝑛p}), for 𝑘∈ {1,…,5}, were obtained also for the 𝐬direction. The objective function 𝛺to be minimized with respect to the vector 𝐩∶= (𝑎1, 𝑏1, 𝑎f, 𝑏f, 𝑎nn, 𝑏nn, 𝑎fs, 𝑏fs)of material parameters in the SEF 𝜓 (Eq. (14)) was formulated as 𝛺(𝐩) ∶= 6 ∑ 𝑗=1 ‖√10 6( 𝝉(𝑗)−𝝉(𝑗)(𝐩))‖2+ 5 ∑ 𝑘=1 (‖ 𝝈(𝑘) f−𝝈(𝑘) f(𝐩)‖2+‖ 𝝈(𝑘) s−𝝈(𝑘) s(𝐩)‖2), (38) where  𝝉(𝑗), 𝝈(𝑘) fand  𝝈(𝑘) sare vectors containing 𝑛pexperimental stress values (e.g.,  𝝉(𝑗)= (𝜏(𝑗) 𝑖∣𝑖∈ {1,…, 𝑛p})), while 𝝉(𝑗)(𝐩),𝝈(𝑘) f(𝐩)and 𝝈(𝑘) s(𝐩)contain the corresponding model values, calculated from (29). We decided to multiply the shear stress errors in (38) by the factor of √10 6(and hence the squares of errors by the factor of 10 6) to eliminate the initial bias arising from the fact that 𝛺(𝐩)is influenced by only 6 simple shear data sets but, in effect, by 10 biaxial data sets (5 for each loading direction). The minimization of 𝛺was realized using the Levenberg–Marquardt algorithm. The same minimization process was applied also to the HO model. The stress relation for the model, analogous to our Eq. (29), can be found in the original paper by Holzapfel and Ogden (2009). The objective function 𝛺HO was defined the same way as 𝛺in (38) and it was minimized with respect to the family 𝐩HO ∶= (𝑎1, 𝑏1, 𝑎f, 𝑏f, 𝑎s, 𝑏s, 𝑎fs, 𝑏fs)of material parameters occurring in the SEF 𝜓HO (Eq. (8)). The Levenberg–Marquardt algorithm always converged toward the same minimizing solution, even from starting points chosen far from the minimum. However, this solution contained a negative value of the parameter 𝑏fs, as can be seen from the first row of Table 1. Negative parameters should generally not be accepted (cf. the end of Section 2) and therefore we decided to modify the defining expression of the objective function 𝛺HO by adding to it a penalty term 8 ∑ 𝑚=1 𝛾𝑚min{0,(𝐩HO)𝑚}2(39) which, by means of the penalty coefficients 𝛾𝑚>0, penalizes any negative component (𝐩HO)𝑚of the vector of material parameters. After this modification, the solution was repeated and we obtained a second minimizing vector (the second row of Table 1) in which, however, 𝑎fs and 𝑏fs were both zero which means that the 𝐼8-term in 𝜓HO was inactivated by the constraints imposed by (39). The stress responses corresponding to both calculated parameter vectors are plotted in Fig. 4against the experimental data. For each curve in the figure, we calculated the coefficient of determination 𝑅2; the obtained values are listed in Table 2together with the mean 𝑅2for each minimizing solution. In addition, the last column of Table 1shows the minima 𝛺HO(𝐩HO)attained by the HO model in both analyses. With our modified version 𝜓 the algorithm did not generate any negative parameters, as can be seen from the calculated values in the third row of Table 1. However, it is also desirable to provide a comparison between the constrained HO model in Table 1and our modified SEF with invariant 𝐾1. For this reason, we decided to perform another solution in which we excluded the term with invariant 𝐼8from 𝜓, thus creating a six-parameter model comparable to the constrained HO model in Table 1. Minimization of the objective function (38) with respect to the reduced parameter vector 𝐩∶= (𝑎1, 𝑏1, 𝑎f, 𝑏f, 𝑎nn, 𝑏nn)then produced material parameters which are listed in the fourth row of Table 1. The stress responses obtained by both fits of 𝜓 are shown in Fig. 5. The 𝑅2values for individual responses are given in Table 2and the weighted errors 𝛺(𝐩)are in the last column of Table 1. It can be seen from Fig. 5that our modified model 𝜓 describes all the experimental sets very well. This is also reflected in the value of the minimum 𝛺(𝐩)in Table 1which is about 80% lower than the corresponding minimum 𝛺HO(𝐩HO)attained by the original model 𝜓HO. The 𝑅2values in Table 2show that the HO model has a good overlap with the simple shear data. It is also possible to achieve an accurate description of biaxial data by assigning sufficiently high weights to biaxial residuals in 𝛺HO, but then the shear responses of 𝜓HO become European Journal of Mechanics / A Solids 111 (2025) 105586 6 J. Vaverka and J. Burša Table 1 Estimated parameters of the strain–energy function 𝜓HO proposed by Holzapfel and Ogden (2009) and of the modified function 𝜓 proposed in this paper. The attained minima of objective functions are denoted as WSSR (weighted sum of squared residuals). Row 1: Unconstrained optimization with 𝜓HO led to 𝑏fs <0. Row 2: Constrained analysis produced 𝑎fs =𝑏fs = 0which, in effect, made 𝜓HO independent of 𝐼8. Row 3: Our model produced a much lower WSSR than 𝜓HO. Row 4: The 𝐼8-dependent term was removed from 𝜓 to make it comparable with the constrained 𝜓HO (row 2); the WSSR is again lower for 𝜓. 𝑎1𝑏1𝑎f𝑏f𝑎s𝑏s𝑎nn 𝑏nn 𝑎fs 𝑏fs WSSR (kPa) (kPa) (kPa) (kPa) (kPa) (kPa2) 𝜓HO 1.287 4.888 2.498 30.266 1.067 28.749 – – 0.615 −14.617 49.094 𝜓HO constrained 1.332 4.760 2.462 30.489 1.031 29.272 – – 0.000 0.000 49.214 𝜓 0.907 7.347 1.298 34.348 – – 0.335 11.111 0.705 0.899 7.520 𝜓 without 𝐼81.074 6.846 1.291 34.962 – – 0.276 11.896 – – 10.124 Table 2 𝑅2values for all curves displayed in Figs. 4and 5. The first five columns, which characterize biaxial tests, contain two values in each cell, where the upper represents the MFD and the lower the CFD. 0.5:1 0.75:1 1:1 1:0.75 1:0.5 (fs) (fn) (sf) (sn) (nf) (ns) Mean 𝜓HO 0.845 0.881 0.872 0.960 0.869 0.988 0.972 0.983 0.977 0.967 0.954 𝟎.𝟖𝟕𝟕 0.583 0.910 0.794 0.774 0.710 𝜓HO constrained 0.843 0.881 0.872 0.959 0.869 0.985 0.972 0.987 0.975 0.968 0.954 𝟎.𝟖𝟕𝟔 0.582 0.909 0.794 0.775 0.715 𝜓 0.924 0.972 0.993 0.990 0.981 0.994 0.984 0.985 0.985 0.993 0.997 𝟎.𝟗𝟕𝟖 0.944 0.988 0.983 0.986 0.952 𝜓 without 𝐼80.922 0.972 0.993 0.989 0.980 0.973 0.985 0.981 0.991 0.976 0.978 𝟎.𝟗𝟕𝟑 0.940 0.985 0.982 0.984 0.944 far worse than those shown in Fig. 4and the material parameters so obtained cannot be used for modeling of shear behavior. The capability of the HO model to describe both biaxial and simple shear data with one set of parameters was tested earlier by Gültekin et al. (2016), Guan et al. (2019) and Martonová et al. (2024) but the obtained fits were also not entirely satisfactory. However, this inaccuracy of the HO model should not be exaggerated because it could be a consequence of the differences in the properties of the biaxial and shear specimens. These differences may be significant despite the fact that Sommer et al. (2015) reduced them by cutting the biaxial and shear specimens from adjacent regions of the heart wall. Guan et al. (2019) and Martonová et al. (2024) also investigated several modified versions of the HO model with improved fitting capabilities but none of them reached the level of accuracy comparable to our model. In particular, Martonová et al. (2024) fitted the same data as we did using modified models with 14, 7 and 5 parameters, and the resulting fits were characterized by the mean 𝑅2of 0.876, 0.894 and 0.872, respectively. Our model outperforms these models because it produces 𝑅2= 0.978 with 8 parameters and 𝑅2= 0.973 with only 6 parameters. 6. Discussion and conclusions The modified model proposed in this paper is, like the original HO model, inspired by the true myocardial architecture as described in several histological studies, some of which were mentioned in Section 2. We wanted to avoid making any excessive changes in the original formulation, thus we preserved its characteristic additively decoupled structure composed of several exponential terms, each of which represents a portion of the constituents of the myocardium. We also maintained the exponential term with coupling invariant 𝐼8, even though its omission or replacement by a simpler term (with only one parameter) might be appropriate, at least for the human data from Sommer et al. (2015) that we used to fit the model parameters. The distinctive characteristic of our approach is that we utilized the anisotropic invariant 𝐾1(whose potential applicability in biomechanics was suggested earlier by other authors) to describe the mechanics of planar myocardial sheets which the former authors described partly through invariant 𝐼4f, and additionally by invariant 𝐼4s. Although our modification consists merely in substituting 𝐼4s by 𝐾1in the SEF 𝜓HO, the results we present in Section 5prove that the effect of such a small change is substantial. It is important to note that this improvement was achieved without increasing the number of material parameters of the model. Unfortunately, there is some controversy regarding the interpretation of the CFD results in terms of the preferred material directions. In some of the previous papers the CFD was identified with the 𝐬 direction (e.g., Holzapfel and Ogden,2009;McEvoy et al.,2018), while in some others with the 𝐧direction (e.g., Gültekin et al.,2016;Guan et al.,2019). From our experience there is usually no in-depth discussion on this issue in the research papers. Usually, the authors seem to arbitrarily choose one of the two options without even mentioning the uncertainty surrounding this topic, which is evident from the inconclusive information and contradictory statements that can be found in the literature. Of course, under such conditions, both viewpoints can be defended by citing a few supportive literature sources. We searched the literature and collected a bulk of relevant information on this matter in Section 4. Based on this information, we found it more appropriate to assign the CFD results to the 𝐬direction. This choice has the additional advantage that it enables to fully exploit the fitting potential of both investigated SEFs, 𝜓HO and 𝜓. If, instead, we associated the CFD with the sheet-normal direction 𝐧, then the terms with 𝐼4s and 𝐾1in 𝜓HO and 𝜓, respectively, would be inactivated in all considered biaxial tests because simultaneous extension in 𝐟and 𝐧directions leads to 𝐼4s <1 and 𝐾1<1. The capability of these terms to improve the quality of fits would thus be considerably limited. One can find the fits of the HO model obtained under such limiting conditions, e.g., in the work of Guan et al. (2019). In order to increase the quality of fits, the authors proposed a modification consisting of a replacement of the invariant 𝐼4s by 𝐼4n. This change, however, makes the biaxial response of the resulting model in the 𝐟 𝐧plane identical to the biaxial response of the original HO model (used by us) in the 𝐟 𝐬plane. The only strain– energy term capable of producing any difference between the models is therefore the 𝐼8-term whose impact, however, is generally very minor, as can be seen from our Figs. 4and 5. This small difference between the two models is perhaps best demonstrated by the fact that the mean 𝑅2of 0.877 that we reached with the original HO model (see Table 2) is almost the same as the value of 0.867 that Martonová et al. (2024) obtained when they fitted the above-mentioned model of Guan et al. (2019) to the same data that we used. Thus, both models are similar, except that the modified version with 𝐼4n exhibits stiffer behavior along the 𝐧direction than along the 𝐬direction, which is unrealistic considering the structural organization of myocardium shown in Fig. 1. If, on the other hand, we admit that the contribution of sheets to the stiffness in the CFD is significant, we can obtain the same results using the original HO model which also correctly assumes that the stiffness European Journal of Mechanics / A Solids 111 (2025) 105586 7 J. Vaverka and J. Burša Fig. 4. Experimental data of Sommer et al. (2015) fitted by the HO model (Eq. (8)). An unconstrained problem was solved first, followed by a constrained problem allowing only non-negative parameter values. The top five panels contain experimental responses to biaxial extension in the mean-fiber direction (MFD) and the cross-fiber direction (CFD), and the corresponding model responses in the fiber direction 𝐟and the sheet direction 𝐬. The title of each graph specifies the applied MFD to CFD strain ratio. The bottom six panels contain experimental and model responses to the six simple shear modes specified in the titles. in tension is highest along 𝐟, medium along 𝐬and lowest along 𝐧. This behavior of the myocardium is also consistent with our model which, however, achieves better fitting results than the HO model with the same number of material parameters. We also want to note that anyone who feels that the CFD should rather be identified with the 𝐧direction (and also does not insist on having 𝐬direction stiffer than 𝐧) can follow the approach of Guan et al. (2019) and simply define 𝐾1in 𝜓 using 𝐬 instead of 𝐧, thereby reinforcing (again) the CFD. The important thing is that one obtains much better biaxial responses with the combination of 𝐼4and 𝐾1than with the combination of two invariants 𝐼4(used by Holzapfel and Ogden,2009;Guan et al.,2019). This underlines the importance of our work because, to the best of our knowledge, the invariant 𝐾1has not been used so far in any constitutive equation for myocardium, which is surprising considering its great potential demonstrated by our results. On the other hand, it must be admitted that there is currently no experimental evidence that the combination European Journal of Mechanics / A Solids 111 (2025) 105586 8 J. Vaverka and J. Burša Fig. 5. Experimental data of Sommer et al. (2015) fitted by the modified model 𝜓 proposed in this paper and by its reduced version obtained by excluding invariant 𝐼8from (14). The top five panels contain experimental responses to biaxial extension in the mean-fiber direction (MFD) and the cross-fiber direction (CFD), and the corresponding model responses in the fiber direction 𝐟and the sheet direction 𝐬. The title of each graph specifies the applied MFD to CFD strain ratio. The bottom six panels contain experimental and model responses to the six simple shear modes specified in the titles. of 𝐼4f and 𝐾1reflects the true structure of ventricular myocardium better than the traditional approach. Another aspect that deserves further discussion is our choice of shear data with a maximum applied shear of 0.5. Regarding the biaxial tests, Sommer et al. (2015) provide complete results only for the maximum applied stretch of 1.1, but they show four different collections of simple shear responses, obtained for the subsequently applied shear levels of 0.2, 0.3, 0.4 and 0.5. The reason for this multistep loading was the pronounced strain-induced softening (Mullins effect) exhibited by the samples in both biaxial and simple shear tests. This means that each increase in load level was associated with an irreversible decrease in stiffness. In other words, different mechanical responses were obtained for different shear levels, which obviously complicates fitting by a hyperelastic model because it is not clear which shear level should be combined with a given set of biaxial data. It would be reasonable to use that shear level in which the constituents most responsible European Journal of Mechanics / A Solids 111 (2025) 105586 9