A Viscoelastic Model for Human Myocardium

Understanding the biomechanics of the heart in health and disease plays an important role in the diagnosis and treatment of heart failure. The use of computational biomechanical models for therapy assessment is paving the way for personalized treatment, and relies on accurate constitutive equations mapping strain to stress. Current state-of-the art constitutive equations account for the nonlinear anisotropic stress-strain response of cardiac muscle using hyperelasticity theory. While providing a solid foundation for understanding the biomechanics of heart tissue, most current laws neglect viscoelastic phenomena observed experimentally. Utilizing experimental data from human myocardium and knowledge of the hierarchical structure of heart muscle, we present a fractional nonlinear anisotropic viscoelastic constitutive model. The model is shown to replicate biaxial stretch, triaxial cyclic shear and triaxial stress relaxation experiments (mean error ~7.65%), showing improvements compared to its hyperelastic (mean error ~25%) counterparts. Model sensitivity, fidelity and parameter uniqueness are demonstrated. The model is also compared to rate-dependent biaxial stretch as well as different modes of biaxial stretch, illustrating extensibility of the model to a range of loading phenomena.


Introduction
The biomechanical function of the human heart is a critical component of cardiac physiology.Beyond the role of its active properties leading to contraction of the myocardium, the passive characteristics of heart muscle play a key role in cardiac pathophysiology, particularly in conditions such as diastolic heart failure [92], heart failure with preserved ejection fraction (HFpEF) [73], and myocardial infarction [14].Patients with these conditions often have poor outcomes due, in part, to patient variability and current challenges in predicting therapy efficacy.Cardiac biomechanical modeling provides a tool for addressing these needs, providing the capacity for patient-specific assessment and model predicted outcomes.These models are playing an increasingly important role in translational cardiac modeling [13] and rely on appropriate constitutive models to predict the passive biomechanical response of the myocardium throughout the cardiac cycle.
Constitutive model characterization of passive myocardial tissue has been a focus of research for over 1.5 centuries [85,7].Experimental studies in animals have shown that myocardial tissue exhibits nonlinear stress-strain response [63], anisotropy in biaxial stretch [18,89,40], and orthotropy under shear [20,21].Recently, these effects were extended and shown in bovine [2,49] and human myocardial tissue [76,77].These experimental insights have driven the development of numerous mechanical constitutive models, with varying degrees of fidelity [13].In most cases, myocardial models have leveraged hyperelasticity theory [34,10,78], defining the stored energy (strain-energy) in response to the loading of muscle tissue.Efforts at developing constitutive relations largely paralleled available data, with early descriptions focusing on transversely isotropic strain-energy equations [39,40,41,25], followed by orthotropic descriptions [16,71,72,36].In the hyperelastic formulation, the transfer from external energy to internal energy (or vice versa) is lossless, providing perfect energetic retention and return.
While current models treat the heart muscle as hyperelastic, this belies the considerable evidence of myocardial viscoelasticity.Viscoelastic response has long been observed in muscle tissue, with early evidence stemming from the work of Blix [7,8,9] who showed hysteresis in ex vivo frog gastrocnemius.This loss was further characterized by Hill and Hartree [31], who demonstrated the loss of energy in the stretch and relaxation of viable/non-viable muscle tissue samples.These viscous elastic effects were initially explained using a spring and dashpot model by Levin [47].Hysteresis and nonlinear stress-strain relations were later demonstrated in the canine papillary by Walker [80] and at the organ-scale in the Langendorff feline heart experiments of Leach and Alexander [44].Viscoelastic relaxation phenomena were observed in the ex vivo beating tortoise heart by O'Brien and Remington [61].Similar experiments were studied in conscious dogs, demonstrating hysteresis and creep in vivo [48].A comprehensive mechanical assessment was later performed in the rabbit papillary muscle by Pinto and Fung [63], showing relaxation, creep, hysteresis along with a modest frequency dependence.Since, viscoelastic behaviors have been reported in many experimental works [18,89,20,21,76,77].
The role of viscoelasticity in myocardial mechanics, while clear experimentally, has yet to be widely adopted the constitutive equations for heart tissue.This is, in part, due to the already complex nature of state-of-the-art cardiac mechanics models -involving multiple nonlinear anisotropic terms with multiple unknown parameters.In addition, the lack of straightforward nonlinear viscoelastic models to build from has limited their extension.A range of rheological analyses have been performed [90], and used as the basis for linear viscoelastic models (see, e.g., [82,24,19]).These efforts were soon extended into nonlinear elasticity theory by Coleman and Noll [15], Truesdell and Noll [60,78], Pipkin and Rogers [64], and extended to quasi-linear viscoelasticity for biological tissues by Fung [24] (see the review by Wineman [83]).Efforts moving nonlinear viscoelastic models into simulations were done by Simo [75] and Holzapfel [33,37,35].Extension of these approaches to the nonlinear viscoelastic behavior in the heart were presented by Huyghe et al. [42] and Cansiz et al. [11].Recently, a study by Gultekin et al. [26] used an analogue to the Maxwell-Wierchert model to characterize viscoelastic effects in the different microstructural orientations, providing the first three-dimensional anisotropic nonlinear viscoelastic constitutive equation for human myocardial tissues.Interestingly, however, the model required significantly different parameter sets and relaxation times depending on the experiment performed by Sommer et al. [76,77].
Here we present a three-dimensional viscoelastic constitutive model framework for the human myocardium.A structural argument is presented based on the hierarchical nature of the myocardial tissue, suggesting the presence of a spectrum of relaxation times.Following the previous works of Simo and others [75,33,37] as well as extensions to fractional approximations [54,32], a fractional anisotropic nonlinear viscoelastic model is proposed that encapsulates phenomena observed in the heart.The model developed for myocardium is fit to recent human myocardial data [76,77], showing mean errors of ∼ 7.65% compared to ∼ 25% for hyperelastic variants.Moreover, the model is used to predict variations of biaxial stretch and stretch rate, showing compelling predictions of the passive muscle response.The developed model uses 11 parameters (compared to 17 − 18 used in [11,26]), which are shown to be unique.This model represents the first nonlinear anisotropic viscoelastic model of human myocardium demonstrated to fit the biomechanical response of myocardial tissue and show rea-sonable predictions of tissue response.
Outlining the work, we begin by reviewing the potential sources of myocardial viscoelasticity, building a structural argument that provides the foundation of the model (section 2.1).Basic notation and hyperelastic formulations for cardiac mechanical models are reviewed (sections 2.2 and 2.3).From these foundations the anisotropic nonlinear viscoelastic power law model for the human myocardium is introduced (sections 2.4-2.6).Sections 2.7 and 2.8 review the human rheological data used in this study as well as the parameterization procedure followed.Results of the model are then presented in section 3, with discussion presented in section 4.

Origins of viscoelasticity in passive myocardium
The complex structure of myocardial tissue has led to discussion over the origins of viscoelasticity.Many constituent components of the myocardium have been implicated as the source of viscoelasticity, including tissue perfusion, extracellular fluid, myocytes, the extracellular matrix (ECM) and others.Simulated poromechanical studies [88] have demonstrated that tissue perfusion can yield stress relaxation.Increased extracellular fluid content is known to significantly influence the biomechanics of tissues [12,27,4].Experimental studies on sarcomeres have demonstrated that the primary contractile proteins of the heart exhibit passive viscoelastic behavior [17].Studies have also shown that the main constituents of the ECM exhibit viscoelasticity [74], suggesting molecular friction as a source of viscoelastic response.While these factors are often considered and advocated for individually, it is highly likely that all factors can contribute to the viscoelastic response of the myocardium with varying degrees of importance depending on the spatiotemporal scales and loading conditions considered.In the following section, we review the evidence for these different factors contributing to the viscoelastic response of myocardial tissue.

Influence of tissue perfusion and extracellular fluid
Due to the surrounding interstitial fluid space and the perfusion of myocardial tissues by the coronary vasculature, it is clear that the heart is a complex poromechanical organ.Debate then arises whether the viscoelastic behavior observed predominantly stems from fluid movement through the porous tissue (i.e. the tissue is nearly poroelastic) or if the solid compartment itself is viscoelastic (i.e. the tissue is poroviscoelastic).Questions also arise regarding the influence in vivo versus the typical testing occurring ex vivo.In tissue studies by [20,21,77,2], non-perfused tissues exhibited significant viscoelastic response.While this could be explained by the movement of interstital fluid, the shear rates required to dissipate the energy observed experimentally would require much larger frequencies and would not explain the seemingly nonlinear loss response observed.This was shown through modeling by Yang [88], who demonstrated that the viscoelastic response due to poroelasticity was not sufficient to explain hysteresis observed in data.However, the presence of extracellular fluid and vasculature has a clear influence on the biomechanics of tissue, with results showing that a change in the aqueous solution directly impacts apparent stiffness [56].Hydration of myocytes and the ECM proteins both have significant impacts on their properties and viscoelasticity, making these factors critical to the passive behavior of tissue.Further, in simulation studies [69], it was shown that pore pressure yields a general stiffening by loading the ECM which could, in turn, influence viscoelastic response of structural proteins.As a consequence, the viscoelastic response of tissue is inextricably linked to the constituents of the extracellular environment whether (or not) the porous flow of fluid plays a leading role in the viscoelastic mechanical response observed.

Influence of cardiomyocytes
The passive mechanics of cardiomyocytes have also been identified as contributors to viscoelastic response.The viscoelastic impact of myocytes is often thought to be significant due to the encapsulated fluid of the cytosol.However, skinned myocytes exhibit strong passive viscoelastic response, thought to arise from molecular friction in the macromolecules such as titin [38].In addition, intracellular structural proteins, such as actin [51], have also exhibited viscoelasticity due to molecular friction.While myocytes exhibit passive viscoelastic response, it is unclear the degree to which passive cellular forces contribute to the total tissue response.Witzenburg et al. [84] demonstrated that decellularization of rat myocardium increased the apparent stiffness of the myocardium 6.7 fold in biaxial tests; an increase that corresponded to the 5.6 fold reduction in cross-sectional area.While this may suggest that cells provide a minor contribution to passive biomechanics, it is clear that the turgid cylinder-like structure of the cardiomyocyte plays a significant role in resistance to shear and compression.

Influence of extracellular matrix structure
The extracellular matrix is acknowledged as a critical, and dominant, component of passive mechanical properties in the heart [81,45,52].The ECM is highly hierarchical [22] (see figure 1), exhibiting important structural organization at multiple scales.A fundamental building block of myocardial ECM are collagen I and III [6] (approximately 300 nm in length).The extremely rigid collagen molecule triple helices are bound to form collagen fibrils, with diameters of 10 − 500 nm and lengths on the order of 10 − 30 µm (see figure 2A).This structural arrangement provides opportunity for significant molecular interaction, yielding potential molecular friction mechanisms as molecules translate relative to one another.Indeed, Shen et al. [74] performed stress relaxation experiments on individual collagen (type I) fibrils, demonstrating viscoelastic response.Characterization of this relaxation response required a 2-component Maxwell-Weichert model, suggesting multiple relaxation times are present at the scale of individual fibrils.It has also been shown that collagen fibrils exhibit nonlinear response, stemming from an uncoiling of end groups leading to progressive recruitment of molecules [23,68].
Collagen fiber formation bundles together many individual fibrils.Fiber bundles vary significantly in diameter (with type I forming larger bundles than type III) and can span over extended distances.To avoid rigid locking between fibrils, a layer of proteoglycans and glycosaminoglycans cover the outer fibril surface [30] (see figure 2B).These large proteins are strongly hydrophilic, binding with water and enabling molecular lubrication between fibrils.Collagen fibers are then linked together to other collagen molecules, often running in close proximity through myocardial tissue (see figure 2C).
At the tissue scale, fibers weave together forming the ECM (see figure 1).Fibers form complex layers of endomysial collagen -surrounding and spanning between individual myocytes -and perimysial collagen -surrounding and spanning between muscle sheets (see figure 2D).Endomysial and perimysial collagen exhibit a complex structure that undergoes predictable molecular realignment and uncoiling under load [66,29,67].Stretch of individual cells requires deformation of the endomysial layer, resulting in molecular friction as fibers translate relative to one another.Scaling up to the level of myocardial sheets, a similar process of deformation and relative translation occurs within the perimysial layer.
Examining the hierarchical structure of the ECM, it is clear that molecular friction and potential drivers of viscoelastic response are pervasive and present at multiple spatial scales.From fibril to fiber to collagen layers, a relative translation of molecules to macromolecular complexes can be observed.The multiscale mechanisms of molecular friction suggest that myocardial tissue is likely characterized by a spectrum of relaxation events occuring at a broad range of time-scales.This is in contrast to some collagen hydrogels which can exhibit less orderly structure [57,58] and can be well characterized by a simple series of relaxation times [50,87,86,43].

Kinematics and notation
Here, we briefly introduce the classic kinematic notation used in nonlinear mechanics (see, e.g., [10,62,34]).Deformation in the heart can be defined by the motion of material points as they move from their reference, Ω 0 ⊂ R 3 , to their physical configurations, Ω t ⊂ R 3 (at some time t ∈ [0, T ]).Marking material points of the reference domain by their coordinate position, X ∈ Ω 0 , the relative motion is defined by the displacement field u : Local deformation of material axes and volumes is characterized by the deformation gradient, F = ∇ X u + I, and its determinant, J = det F > 0, respectively [34].Often, the heart is approximated as an incompressible tissue (e.g., J = 1), though debate about this exists in the literature [3].Measures of stretch are given by the right and left Cauchy Green tensors, defined as C = F T F and B = FF T , respectively.
Often it is convenient to define constitutive relations using invariants of stretch / strain tensors (here generically denoted by A) given by [10] It can also be convenient to consider isochoric invariants, replacing A with Ā = A/III 1/3 A in Eq. ( 1) (e.g., ĪA = Ā : I), as these can facilitate the isochoric / deviatoric split in hyperelastic formulations.
Considering anisotropic materials, such as myocardium, a common approach is to define mutually orthogonal local microstructural directions, e.g., {e f , e s , e n }, at each material point.Here e a denotes the unit orientation vector along myofibers (f), sheets (s) and sheet normals (n) [46].The physical orientation of these microstructural directions is given by ε a = Fe a , a ∈ {f,s,n}.Hence, the squared stretch along microstructural directions can be expressed by the pseudo invariants, Generalizing Eq. ( 2), where sym( A) = 1 2 (A + A T ) is the symmetric transformation.These invariants play an important role in defining the constitutive behavior in myocardial tissue [36].

Hyperelastic formulations of the myocardium
Though a number of constitutive equations have been developed to describe myocardial tissue [71,72,36], in this study we focus on two of the most common, the Holzapfel and Ogden [36] (HO) and the Costa [16] models.The HO model is a hyperelastic incompressible strain-energy function, widely employed in modeling studies and is appealing for both theoretical reasons (e.g., convexity and objectivity) as well as parameter identifiability [72,28].Recalling briefly, the incompressible HO strain-energy W e : R 3×3 × R → R ≥0 can be written as where and δ kl denotes the Kronecker delta (zero unless k = l, in which case its 1).In this formulation, the strain-energy is written as a sum of components.The isotropic term, W iso , provides the isotropic bulk distortional energy associated with tissue deformation.Anisotropic terms W kl (where k,l ∈ {f,s,n}) are introduced to account for the varying distortional energy associated with deformation along microstructural directions.We note that the anisotropic invariants I ff and I ss are thought to not support compression, with W ff and W ss set to zero when I ff , I ss < 1 [36].
In the HO model, the set of anisotropic terms depend on 8 fitting parameters (denoted by a's and b's).In the HO model, the second Piola Kirchhoff stress (PK2) S can be written as with S = {ff, ss, fs} and S p = pJC −1 .Through push forward operations (i.e., σ = J −1 FSF T ), the second Piola Kirchhoff stress tensor can be extended into the Cauchy stress tensor.We note that the original HO formulation considered the strain-energy defined using standard invariants [36], with later forms defined using isochoric invariants and/or dispersion [26].
Another constitutive equation that is widely employed is the Costa model [16].This model is well-posed in terms of convexity and objectivity [36] and relies on an orthotropic formulation of the exponential Fung-type law.While typically written in terms of rotated Green Lagrange strain tensors, the Costa model can similarly be expressed in terms of invariants of the right Cauchy Green strain tensor, e.g. where and S = {ff, ss, nn, fs, fn, sn}.This form is useful for comparison of these models, illustrating that the Costa model shares the same underlying exponential invariant formulation; however, strain-energy terms are grouped in one term W(C), through multiplication.In this case, the Costa model is comprised of 7 material parameters, an outer scaling constant, C, along with anisotropic scaling constants, b.The PK2 stress tensor can be written as Both HO and Costa models utilize exponential forms commonly employed in collagenous biological tissues [23,68] reflecting the toe, heel and linear regime of collagen fibrils and the progressive recruitment of coiled fibers with stretch.Further, while in this paper we focus on incompressible forms, both models have been used considering the myocardium as compressible (or nearly incompressible), where S p is replaced with an appropriate compressibility term.The primary variation between models is the separable form employed in HO compared to the coupled exponential form of the Costa model.Both model forms have advantages and disadvantages in terms of model performance and parameter identifiability.

Extension to cardiac viscoelasticity
Viscoelastic phenomena in muscles have been described using a range of linear spring, dashpot and spring pot system models [47,55].One of the simplest models that captures stress relaxation (impulse stretch) and creep (impulse stress) is the Zener model, containing a Maxwell model in parallel with a spring (cf.[34,19]).Denoting the stress in the elastic branch of the model as σ e 0 (t) = E 0 ε(t), the total stress in the Zener model can be written as where τ = (E/η) −1 is the relaxation time and B = E/E 0 is the relaxation weight [90].In this model, the asymptotic elastic response is given by the first term, while the viscoelasticity is encapsulated in the second term.Eq. ( 11) provides a template for defining the stress response [1].Extending this idea to nonlinear materials, a straightforward approach is to replace the elastic stress component σ e 0 with its nonlinear hyperelastic variants.Replacing σ e 0 with an appropriate nonlinear PK2 form S e , yields the quasi-linear viscoelasticity (QLV) Fung model, where, for the Zener model, K(s − t) = B exp{(s − t)/τ}.Here the QLV Fung model extension is considered linear due to the linear dependence on S e .As noted in section 2.1, myocardial tissue is strongly influenced by viscoelastic phenomena across a hierarchy of scales, with collagen fibrils themselves exhibiting multiple relaxation times [74].A straightforward adaptation is to consider a nonlinear variant of the Maxwell-Wiechert model (or Generalized Maxwell model) [79,70].Considering n viscoelastic elements in parallel, due to linearity, this could be modeled simply by altering the QLV form, with the relaxation Here, B k and τ k represent the relaxation weights and relaxation times of the various viscoelastic elements.Note that, here, we have assumed that the branches share the same nonlinear form S e , which need not be the case.While this form provides generality for encapsulating the viscoelastic response of a material in the QLV framework, it also requires unique determination of all the relaxation weights and times, which is challenging to obtain experimentally particularly for nonlinear anisotropic materials like the myocardium.However, if these weights and time constants were distributed in a well-characterized way, it's possible that the spectrum could follow a simplified approximate form.
In myocardial tissue, viscoelastic response has been observed stemming from molecular friction processes, flow or intracellular effects (see section 2.1), all of which occur at one or more spatiotemporal scales.In the context of myocardial ECM, friction processes are observed at the basic scale of the fibril and extend, hierarchically, as we consider different scales within the ECM (see figure 2).To argue for a specific form, we consider a representative volume of myocardial tissue.At the microscale, the relaxation times are relatively small reflecting the fact that friction occurs between molecules at small spatial scales.Bridging toward larger length scales, the relaxation times increase reflecting that friction occurs between larger conglomerates of molecules or whole tissue structures.The effective density of these different mechanisms also varies across scales.Within a representative volume, the occurrence of microscale phenomena is abundant, while the occurrence decreases toward larger spatial scales.Under these heuristics arguments, the composite relaxation weight could be envisioned as a power law distribution, i.e.K(z) ∝ z −α (see figure 2).

Extension to fractional viscoelasticity
Based on the argument that, due to the hierarichal structure of myocardial tissue, the relaxation function K resembles a power distribution, the QLV model simplifes to the nonlinear fractional Kelvin-Voigt model, i.e.
where D α t denotes the α-order Caputo derivative [65] and B 0 is a strictly positive constant.Here, B 0 modulates the relative strength of the viscoelastic response, while α modulates the distribution of relaxation times.Indeed, for α = 0 the relaxation spectrum becomes infinitely dependent on instantaneous time-scales, reducing the term to a pure elastic response.In contrast, as α tends to 1, then the fractional operator becomes a nonlinear dashpot.Use of this model significantly reduces the complexity and number of parameters required to construct the material response.This simplification is not for free; however, but instead reflects an inbuilt assumption that the relaxation spectrum exhibits power-law behavior.Generalizing this framework, the viscoelastic PK2 stress tensor can be conveniently written as a fractional viscoelastic differential equation, e.g.
or the analogous form, where S v and S e denote the nonlinear viscoelastic and elastic responses, and δ > 0 is a scaling weighting the history dependence on the PK2 stress tensor itself.Note, this is a departure from the linearity seen in the QLV model.Eq. ( 16) was shown to exhibit realistic viscoelastic behavior across a range of testing scenarios [90], providing advantages to the fractional form in Eq. ( 14).Eq. ( 17) performs similar to Eq. ( 16), with the primary difference that the model retains a purely elastic term S e .

A fractional viscoelastic model for the myocardium
In this paper, we propose a viscoelastic model capable of capturing the behavior of myocardial tissue.Building from sections 2.3 and 2.4, we propose the model form, In this case, the elastic form S e denotes the response of the underlying base structure of the tissue and is characterized by a simple neo-Hookean model with a single stiffness parameter, a, e.g.
The complex nonlinear and viscoelastic response of the tissue is incorporated through S , which is governed by the fractional differential equation In this formulation, the underlying material response is dictated by the choice of S v , with the viscoelastic effects controlled by the parameters α and δ.Mimicking the structural forms of both HO and Costa models, we choose S v to be, where S = {ff, ss, nn} and S ⊥ = {fs, fn, sn} and Like the HO model, the first term contains an isotropic scaling, W 1 ; however, this term scales weighted structurally anisotropic terms.As a hyperelastic contribution, this term would lead to a diagonal matrix at the reference state, suggesting the reference state is not a stable state.However, as part of the fractional viscoelastic model (differentiated by the Caputo derivative), this is no longer a factor.Similar to the Costa model, the invariant terms are grouped, with the first term covering the sum of the diagonal invariants (note, I C = I ff + I ss + I nn ) and the second incorporating the shear invariants.Similar to the first term, W 2 scales weighted shear terms.Unlike the Costa model, weighting is done through scalings b kl which are not reflected in the exponents of W 1 or W 2 .While incorporation of scalings b kl into the exponential scalings of S v would make it straightforward to write as a strain-energy function, given the model is viscoelastic, this constraint is lifted.Moreover, grouping terms with two exponents, b 1 and b 2 , significantly simplifies the parameterization process.Last, the model considered in this paper utilizes an incompressible formulation, where S p = pJC −1 .However, other formalisms are possible to incorporate compressibility (though care must be taken when considering these forms [59]).Important properties for constitutive models are objectivity and material frame indifference [83].While these properties are well-known for S e and S p , these properties must hold for S .Examining Eq. ( 20), it's sufficient to require that these properties hold for S v .Specifically, that S v is invariant to rigid (potentially dynamic) rotations of the physical frame and that, given an orthonormal rotation H of the material frame, HS v (C, e m ⊗ e n , . ..)HT = S v (HCH T , He m ⊗ e n H T , . ..). (22) Objectivity holds as S V is a strict function of the invariants of C. Further, material frame rotations leave all invariants unchanged

Experimental testing of human myocardium
This modeling study relies on rheological tests performed on human myocardial samples by Sommer et al. [77,76].While complete details on the experimental data acquisition can be found in [77,76], we reiterate the basics of data acquisition here (see figure 3).Briefly, human heart muscle samples were collected during transplant surgery, infused with 200mL cardioplegic solution (CPS, Celsior by Genzyme Corporation), inserted into a path of 1000 mL of CPS and cooled to 4 • C. From the collected hearts, samples were cut into 25 × 25 mm 2 thin sheets for biaxial testing and (4 mm) 3 cubed samples for shear testing.Passive tests of the cardiac muscle were performed at 37 • C, in a bath of CPS with 20 mM 2,3-butanedione monoxime (BDM) .
Triaxial shear tests were performed on cubes extracted from adjacent locations of the biaxial samples.The cube sides were cut in order to align with the tissue microstructure (figure 3a).Shear deformations between 10 and 50% were performed, with a 10% increment and loading speed of 1 mm/min.Shear relaxation tests were also performed at 50% at a ramp speed of 100 mm/min and a duration of 5 min.As a cube can be used to test two orthogonal directions, a total of three specimens were needed to test all six modes of simple shear.For triaxial testing, two deformation cycles were used for preconditioning, with the data being acquired during the third cycle.

Parameter identification and analysis
The proposed model parameters were fit to human triaxial and biaxial experimental data collected in [77], assuming idealized shear and biaxial kinematics, i.e.F = γe k ⊗ e l + I, and where e k , e l denote shear along different microstructural directions, and γ, γ 1 , γ 2 denote the time dependent amount of stretch or shear applied for each test (see figure 3).Here we refer to the various data tests as groups -relaxation, cyclic shear and biaxial stretch, which will be indicated by superscripts i ∈ {r, c, b}, respectively.The relaxation and cyclic shear groups comprise six tests each, indicated by subscript kl corresponding to the shear directions kl ∈ M s = {fs, fn, sf, sn, nf, ns}, while the biaxial stretch group comprises 2 tests indicated by subscript kl ∈ M b = {ff, ss}.Hence, a total of 14 tests were used for model fitting.For each of those, let σ i kl denote the respective component of the Cauchy stress, with σkl marking the values measured in the experiments, and σ the values computed using the model.The bold versions, σ i kl and σi kl , denote the stress over time for the i th test set and kl th stress component.Note, all models were preconditioned following the data protocol, and compared at their final cycles (where applicable).

Minimisation problem
In this paper, we consider the fit of three models: the HO model, Costa model, and proposed viscoelastic model from section 2.6.All models considered were fit to data by minimizing the objective function across all 14 tests simultaneously (Eq.( 23)) to determine the model's set of N-parameters (denoted as θ) θ = arg min where the objective function, χ, is given as and R i kl = σi kl 2 is introduced to give tests equivalent weights (irrespective of the magnitude or number of time points).Note that χ gives the relative error across all tests, with χ = 1 denoting 100% error.In this minimization, a parameter set is found to best match all datasets.
In this study, the data is representative of the myocardium behavior in relaxation, cyclic shear and biaxial stretch, but it is not guaranteed to come from the same sample (e.g., for the same shear mode in relaxation and cyclic oscillations), or even from the same heart.Therefore, innate variability is introduced that may be reflective of animal to animal variability.Here we assume that the relative magnitudes test to test may vary as a result of different samples, but the shapes should be maintained.This is achieved through the introduced relative scalings β, which are selected to minimize the difference between model and data across all time points.In order to meaningfully preserve the anisotropy, β = (β r , β c , β b ) was applied across all kl directions for given tests.The selection of β was done iteratively with an initial guess β = 1, solving the minimization for y, minimizing for β and repeating until convergence.
Errors were also quantified on a test by test basis to understand areas of model strength and weakness using the related test-specific objective function, where the β i and θ are determined from the minimization in Eq. (23).Minimization for the Costa and HO models were done using the nonlinear lsqnonlin minimization routine in MATLAB.Multiple initial guesses were selected for the model parameters, with the reported results selected using the best fit.
To examine the behavior and uniqueness of the parameter space of the proposed viscoelastic model, a parameter sweep was performed.To generate model predictions, each experiment had to be simulated based on the given kinematics to solve the fractional differential equation forward in time (see [91,90] for details).In this case, we exploited the linearity of some model parameters θ l , which could be found through a simple least squares matrix solve for a given the set of nonlinear parameters θ * ; reducing the parameter space from R 11 + to R 4 + .Uniqueness of the linear parameters was then determined by the solvability of the least squares system.Because of the simple kinematics and known stresses during the relaxation tests, it was also possible to reduce the parameter space to R 3 + by determining the functional relationship between δ and α during the relaxation tests.Sweeps were performed over b 1 ∈ {0, 1, 2 . . .20}, b 2 ∈ {0, 0.5, 1 . . .10}, and α ∈ {0.2, 0.22, 0.24 . . .0.4} resulting in 4, 851 minimizations and 67, 914 test simulations.Note, a coarser sweep over higher parameter values showed continued growth in χ.This sweep indicated that the minimum is achieved at α = 0.24, b 1 = 9 and b 2 = 2.A refined sweep was then conducted around these values, as follows: α ∈ {0.2, 0.21, 0.22 . . .0.28}, b 1 ∈ {7, 7.5, 8 . . .11} and b 2 ∈ {1, 1.25, 1.5 . . .3.5}.

Model sensitivity analysis
Sensitivity was assessed from two perspectives: the effect of noise in the data on the best-fit material parameters, and the effect of parameter perturbations on the fitting error.In the first instance, 100 noisy datasets were produced via adding unbiased uniformly distributed noise vector over time η, to the measured stress values: σi kl = σi kl + η.The noise level was set to 10% of the peak stress across a test for a specific mode, e.g., For each noisy dataset σi kl , a model fit with a different set of optimal parameters θ were obtained, following the same parameter fit procedure.Next, the optimal parameters θ min were perturbed by ±10% of the original value, by turn.The error function χ was computed for each local perturbation, according to Eq. (24).

Results
Figure 4 shows the best fit for HO, Costa and proposed viscoelastic (VE) models to the data.Example curves for the relaxation, cyclic shear and biaxial tests are shown along with bar plots showing test-specific errors for each of the 14 datasets, computed according to Eq. ( 25).The overall errors, computed using Eq. ( 24), were 24.17, 24.95 and 7.65% for HO, Costa and VE models, respectively.In addition, figure 5 explores the VE model error across the α (0.2 to 0.4), b 1 (0 to 20) and b 2 (0 to 10) parameter space.The error across all datasets, computed according to Eq. ( 24), is shown in the bottom right panel.The error computed separately for the relaxation, cyclic shear and biaxial tests is shown in the top left, top right and bottom left panels, respectively.The isosurfaces help visualizing the error behavior in the 3D space.In each of the four panels, the minimum error is indicated by the white sphere with the overall minimum uniquely identified at α = 0.23, b 1 = 9.5 and b 2 = 2.25 over all datasets.The errors over the individual tests were computed using Eq.(25).  Figure 7 shows the prediction of the models under biaxial stretches of various ff:ss ratios (1:1 -blue, 1:0.75 -cyan, 0.75:1 -green, 1:0.5 -red and 0.5:1 -black).The models are shown in the top left (HO), top right (Costa) and bottom left (VE) panels, while the original data, from [77], are presented in the bottom right panel.The solid curves indicate the fiber response, and the dashed curves show the response in sheet direction.Note that here none of the models was fit to data (except for ff:ss ratio of 1:1 which was part of the original fitting dataset), and the curves represent predictions obtained by using the parameters shown in figure 4.
Figure 8 shows the prediction of the viscoelastic model in ff:ss biaxial stretch (left), at three frequencies: 0.01, 0.033 and 0.1 Hz.The corresponding data, from the original study [77], is shown in the right panel.At 0.01 Hz, the testing conditions are identical to the biaxial test used for the model fitting.However, the data used for the fitting and the data shown in figure 8 come from different samples, with the data in the frequency test reaching fiber and sheet peak amplitude smaller by a factor ∼1.48.
Figure 9 shows the VE model prediction to shear compared with data from the original study [77] on human myocardium as well as porcine myocardium.The shear tests conducted for human myocardium start from low (γ = 0.1) to high shear (γ = 0.4), while the data from porcine myocardium start from high shear (γ = 0.4) to low (γ = 0.1).In both cases, 4-6 cycles are performed at a given shear level with the final cycle shown.Predictions from the VE model are shown (center) using the     shear at multiple levels from human myocardium [77], (center) model predictions of shear response, and (right) cyclic shear in porcine myocardium.Experiments used the same rig and protocols, with the only exception that (left) increased shear levels (preconditioning each level), while (right) started from higher shear, decreasing shear levels.

Analysis of the proposed viscoelastic model
The model proposed in this work introduces viscoelasticity through a fractional approach.This choice reflects a hierarchy of relaxation times mirroring the hierarchical structure of myocardial tissue and the different spatiotemporal scales that lead to viscous dissipation.The proposed viscoelastic model presented here relies on an underlying form that combines aspects of the HO and Costa hyperelastic models.As in the HO model, the model terms used here are based on invariants that reflect the microstructural composition of the passive myocardium.However, the terms are not completely separated, allowing for coupling between invariants, as seen in Eq. (21).Coupling of the microstructural directions is characteristic of the Costa model, where all stretches are inherently coupled.Overall, the VE model captures the characteristics of the data, as seen in representative examples of the three deformation protocols in figure 4. In relaxation, the model can capture the initial peak and the subsequent decay.In cyclic shear and biaxial stretch, the hysteresis and nonlinearity are well matched.
The viscoelastic model entails 11 parameters uniquely determined through the fitting process.Four of these have a nonlinear effect on the model -the fractional order α, exponential powers b 1 and b 2 and the viscoelastic PK2 scaling δ, while the remaining 7 parameters act as linear scalings of terms.Parameter δ is determined from the relaxation tests, at each α value, as explained in section 2.8.For the other 3 nonlinear parameters, a sweep is carried, and at each point in the 3D space a set of 7 unique linear parameters can be identified.Figure 5 shows the error behavior across the 3D space of the sweep.Overall, it can be observed that the region of minima for the cyclic shear and biaxial tests (top right, bottom left) is close to that of the entire dataset error (bottom right), with α between 0.23 and 0.28, intermediate b 1 and small b 2 .For the relaxation tests, the α value is similar, yet b 1 is small and b 2 is large.However, of all groups, the best parameter region in the relaxation tests appears the shallowest, with errors being small as long as α < 0.25.This indicates that the fit of the exponential powers b 1 and b 2 is driven by the cyclic shear and biaxial tests.Among the three groups -relaxation, cyclic shear and biaxial, the smallest errors are achieved in relaxation, as it can be seen in both Figs. 5  and 4.
Sensitivity analysis was conducted on the viscoelastic model and its parameters.Compared to the original values (figure 4), the noise led to at most 4% standard deviation and averages close to the reference values, as seen in figure 6.This suggests that the model parameters are identifiable and robust to noise.A perturbation of ±10% in b ff and b ss leads to the model error increase of ∼ 1.5%, while ±10% perturbation of α and b 1 leads to a significant error increase -3 to 11%.Importantly, all parameters alter the objective function, showing the model parameters are observable.

Model-based predictive response
To examine the response of the proposed VE model, predictions of other tests not used for training were simulated, including biaxial stretch to different fiber/cross-fiber ratios (figure 7), frequency response (figure 8), and different degrees of maximal cyclic shear (figure 9).Applying biaxial stretch at different ratios, the proposed VE model shows excellent predictive behavior by capturing the inherent nonlinearity and hysteresis.Importantly, it also demonstrates the inherent coupling of fiber/cross-fiber stretches seen in the data, whereby change in the stretch applied in one direction influences the stress in other directions.Another predictive behavior of the VE model is shown in figure 8, for biaxial stretch tests at 0.01, 0.033 and 0.1 Hz.The data show that changes in frequency of one order of magnitude yield a modest increase in hysteresis and peak stress amplitude (up to ∼ 19%).Modest increase in the viscoelastic model is also observed (∼ 41%), due to the fractional differential form.While the predictive increase is higher, this model captures this modest growth and would likely perform better with multi-frequency data created at higher frequencies.
Lastly, the VE model behavior is investigated at lower cyclic shear levels.From the original study data [77], strain softening can be observed as the shear increases with fairly sustained hysteresis, see figure 9 (left plot).In contrast, the proposed VE model (center plot), predicts decreased hysteresis at lower shears and nested curves that do not exhibit strain softening.However, the original data was acquired following a protocol of progressive increases in shear strain, with preconditioning cycles at each level.This means that the sample was not preconditioned to the largest shear levels until later cycles.To demonstrate this impact, the test was repeated using porcine myocardium (using the same experimental protocol) whereby the shear protocol was stepped from largest to smallest shear strain (right plot).In this case, the passive tissue response resembles that predicted by the VE model.

Comparison with other models
To investigate the relative impact of the proposed VE model, results were compared with standard hyperelastic HO and Costa models (see figure 4).As expected, both hyperelastic models exhibit a constant stress in relaxation and no hysteresis; however, both models do well at capturing the behavior of the data with errors of 24.17 and 24.95%.It can be seen that the HO model performs better in describing the nonlinear trend of the biaxial test, with the errors in this group being ∼ 10% smaller than for the Costa model.In contrast, the proposed model reduces the error metric by approximately a factor of 3 to 7.64%, with the largest reduction in relaxation tests.While these results are expected to improve (in part, due to the increase in parameters from 7-8 to 11), these tests provide an important benchmark for understanding the potential benefit of using a viscoelastic modeling approach.
Another important comparison is with the previous viscoelastic model published by Gültekin et al. [26] that fits the data utilized here.This paper extended the HO model using a nonlinear Maxwell approach [90] that relied on 18 parameters fit to each test separately.Because this model utilized a single relaxation time, the best fit to the relaxation data showed challenges capturing the entire decay spectrum.In contrast, the proposed VE model captures this decay spectrum while also fitting other datasets.The model also shows difficulty reproducing the cyclic shear data, particularly compared to the model proposed.In addition, the model fitting for the Gültekin model requires determination of 14 nonlinear parameters, significantly increasing the cost of parameterization and analysis.
Examining the predictive responses, a key observation comes when comparing the biaxial data at varying levels of stretch, see figure 7.In this figure, the classic HO model shows an uncoupling of stretch along microstructural directions, whereby loading the fiber to the same stretch and varying the stretch across fibers produces nearly the same load.This is in contrast to the Costa model and the experimental data, which show that these loads are not fully independent.In contrast, the VE model provides varying responses that qualitatively match the behavior of the data.In our testing, adding fiber dispersion to the HO model does show improvement; however, this tends to come at the expense of accuracy in other tests and adding these parameters did not substantially improve the model response.As this model forms the basis of that presented in [26], it is likely that similar challenges to biaxial prediction would be observed.

Study limitations
In this study, data from 14 tests were utilized to parameterize all models -a challenge which is often not attempted within constitutive model studies.While integration provides arguably a more complete result, the challenge comes that no tests stem from the same samples.Hence, the inherent variability in tissues make the analysis challenging.To circumvent this, all models were arbitrarily scaled based on testing groups (relaxation, cyclic shear, biaxial stretch).This allowed models to capture the essential shape, without being heavily biased by the total amplitude.An improvement to this study could be the use of new testing rigs [2], capable of providing a wide range of tests on a single sample.

Conclusions
In this paper, we present a nonlinear viscoelastic constitutive model for passive myocardium using a fractional approach.The model is tested against human myocardial data across a range of testing protocols, demonstrating effective reproduction of experimental measurements as well as strong prediction of other tissue measures.The model builds on the idea of cardiac viscoelasticity stemming from its hierarchical structure, yielding a spectrum of viscoelastic phenomena that space spatiotemporal scales.The model is compared with classic hyperelastic models, demonstrating a significant improvement in fitting (mean error ∼ 7.65% compared to ∼ 25%), as well as in predictive response (particularly for variable biaxial stretch).The model also is shown to exhibit a unique set of parameters that are observable and it is robust.The proposed VE model presents one of the first constitutive models aimed at capturing viscoelastic nonlinear response across multiple testing regimes, providing a platform for better understanding the biomechanics of myocardial tissue in health and disease.

Declaration of competing interest
The authors declare no competing interests.

Figure 1 :
Figure 1: Scanning electron microscopy of ECM structure in the (A) rabbit and (B,C) canine myocardium.(A) Illustrates the detailed microstructure and collagen fibers composing the endomysial collagen layer that normally surrounds and interlinks cardiomyocytes.(B) Shows magnification of the ECM structure, showing pockets normally encasing multiple myocytes as well as the coronary microvasculature.(C) Magnification of the ECM, showing myocardial sheets separated by perimysial collagen layers (marked with arrows).Reproduced with permission from [53, 5].

Figure 2 :
Figure 2: Schematic representation of the hierarchical structure of the extracellular matrix, and its representation as a power spectrum.(A) Illustration of the collagen triple helix and structural arrangement of the collagen fibril.(B) Illustration of the lattice arrangement of fibrils within a collagen fiber, showing fibrils, proteoglycans and their hydration.(C) Illustration of bundled fibrils within collagen fibers and their cross section.(D) Illustration of the myocardial tissue showing myocytes and capillaries surrounded by endomysial collagen fibers and encased in sheets covered by perimysial collagen.(Bottom Row) Illustration of multiscale friction processes yielding a fractional relaxation response moduli.(Left) Sketch of multiscale friction processes (A) within fibrils, (B) between fibrils, (C) at the endomysial collagen scale, and (D) at the perimysial collagen scale.(Middle) Effective density of relaxation phenomena within a representative volume.(Right) Combined multiscale relaxation response modulus with colors indicating different scales of response and its representation using a fractional model (red curve).

Figure 3 :
Figure 3: Illustration of the experimental tests performed on human myocardial samples in [77].(a) Reference tissue block with shear deformations applied in different directions relative to the underlying tissue microstructure.(b) Reference biaxial sample with stretch applied in fiber and sheet directions.(c) Experimental test rigs for triaxial and biaxial experiments.

Figure 6
Figure 6 illustrates synthetic noisy datasets alongside the original data, showing the resultant noisy data for cyclic shear (fs, top left) and biaxial (top right) stretch (ff).Equivalent datasets were generated for all test cases, and replicated 100 times to examine the sensitivity of the fit to noise.For each of the 100 noisy datasets, a new set of VE linear parameters θ l were determined.The parameters' mean (red marker) and standard deviation (error bars) are shown relative to the original parameters, indicated by the baseline 1 (bottom right).To further understand model sensitivity to parameters, variation of the objective function was found when changing parameters by ±10% (bottom left) for both linear and nonlinear parameter sets.The original parameter values can be seen in figure 4. The dotted line indicates the minimum model fit error of 7.65%.

Figure 5 :
Figure 5: Values of the objective function χ across fractional order values α ∈ [0, 1] and exponential powers b 1 , b 2 ∈ [0, 20].The heat maps indicate how the error changes across the 3D space.Isosurfaces are aiding in the visualisation of the error, with the white spheres indicating the location where the minimum is achieved.(Top Left) Error χ r across the relaxation tests only.The minimum (2.86%) is achieved at α = 0.2, b 1 = 1 and b 2 = 7. (Top Right) Error χ c across the cyclic shear tests only.The minimum (9.65%) is achieved at α = 0.28, b 1 = 8.5 and b 2 = 1.(Bottom Left) Error χ b across the biaxial stretch tests only.The minimum (7.90%) is achieved at α = 0.26, b 1 = 10 and b 2 = 0. (Bottom Right) Error χ across all tests.The minimum (7.65%) is achieved at α = 0.23, b 1 = 9.5 and b 2 = 2.25.

Figure 6 :
Figure 6: Sample datasets for the fs and ff modes of deformation with added 10% unbiased noise compared to the original experimental data.(Top Left) Noisy data in fs cyclic shear and (Top Right) ff biaxial stretch.(Bottom Left) Relative variation of model parameters obtained by fitting 100 noisy datasets.For each parameter the mean (shown as the red dot) and the range of one standard deviation (solid black line) can be compared relative to the original value (dotted line).(Bottom Right) Variation of the objective function χ due to 10% variation in parameters.

Figure 8 :
Figure 8: The behavior of the proposed VE (left) against data from the original study [77] (right), in biaxial stretch tests with variable loading frequency (0.01 Hz -green, 0.033 Hz -red, 0.1 Hz -blue).Solid curves show the behavior in the fiber direction, and dashed curves in the sheet direction.The model employs the parameters shown in figure 4 and β b = 0.45.

Figure 9 :
Figure 9: Qualitative comparison between data and model in cyclic shear.(Left)shear at multiple levels from human myocardium[77], (center) model predictions of shear response, and (right) cyclic shear in porcine myocardium.Experiments used the same rig and protocols, with the only exception that (left) increased shear levels (preconditioning each level), while (right) started from higher shear, decreasing shear levels.
DN acknowledges funding from the Engineering and Physical Sciences Research Council (EP/N011554/1) and the Engineering and Physical Sciences Research Council Healthcare Technology Challenge Award (EP/R003866/1), and support from the Wellcome Trust EPSRC Centre of Excellence in Medical Engineering (WT 088641/Z/09/Z) and the NIHR Biomedical Research Centre at Guy's and St.Thomas' NHS Foundation Trust and KCL.The views expressed are those of the authors and not necessarily those of the NHS, the NIHR, or the DoH.