A multiphysical computational model of myocardial growth adopted to human pathological ventricular remodelling

We present a novel three-dimensional constitutive model that describes an electro-visco-elastic-growth response on the myocardium with a fully implicit staggered solution procedure for the strong electromechanical coupling. The novel formulations of the myocardium allows us to simulate and analyze the remodelling of actively contracting human ventricular heart models which consist of growing viscoelastic myocardium where the growth direction is determined based on its mechanical state at each time step. The total deformation gradient is multiplicatively decomposed into a mechanical-active part and a growth part, where the mechanical-active part is further split into elastic, viscous, and active components. Unconditional stability of time integration is ensured by a backward Euler integration scheme. With the developed model, the myocardium can experience stretch-driven longitudinal (fibre) growth and stress-driven transverse (cross-fibre) growth. To validate the developed approach, two simulations regarding pathological ventricular remodelling are implemented: two divergent types of remodelling of a left ventricular model driven by hemodynamic overloads and ventricular remodelling triggered by acute myocardial ischemia in a biventricular heart model.


Introduction
Left ventricular (LV) remodelling is the process of structural and functional changes in response to cardiovascular damage and a precursor of clinical heart failure [1].In principle, LV remodelling can be largely categorized as physiological and pathological.The former is an adaptive remodelling (healthy) and encountered, for instance, in athletes or pregnants, which results in improved cardiac function and is reversible when the subject is no longer exposed to the elevated hemodynamic load [2].Pathological LV remodelling is, however, a maladaptive remodelling process (unhealthy) that harms the cardiac function, mostly irreversible and eventually might B Michael Kaliske michael.kaliske@tu-dresden.de Yongjae Lee yongjae.lee@tu-dresden.deBarış Cansız baris.cansiz@tu-dresden.de 1 Institute for Structural Analysis, Technische Universität Dresden, 01187 Dresden, Saxony, Germany trigger heart failure and life threatening cardiac arrhythmia [3].
A prevalent pathology of remodelling can be developed after (acute) myocardial ischemia [4,5].In most cases, myocardial necrosis and disproportionate wall thinning are observed in the region of ventricular ischemia [6].The weakened region is unlikely to carry the pressure and volume load applied on the ventricles as the other healthy myocardium outside of the ischemic area, which causes consequently local dilatation in the corresponding tissue.In case that the cardiac ischemia vanishes with cells recovered in its initial stage of development, the cardiac function is maintained.However, if the heart undergoes the ischemia in the long run, it is exposed to the risk of pathological LV remodelling that progressively leads to decline in cardiac function.Hemodynamic overload is also common cause of pathological ventricular remodelling, i.e., pressure and volume overload, which leads to different patterns of ventricular remodelling.Pressure overload causes a myocardial growth due to additional sarcomeres in parallel so that ventricular wall thickening is observed at organ level, which is denoted as concentric remodelling [7,8].On the other hand, volume overload leads to growth in cardiomyocyte length through additional sarcomeres in series.This type of myocardial growth induces LV dilatation at organ level, which is denoted as eccentric remodelling [7].
The theory of finite growth can be traced back to Rodriguez et al. [9], who introduced an incompatible growth configuration that multiplicatively decomposes the deformation gradient into an elastic and a growth part.Actually, it was developed for finite strain plasticity [10] and constitutes the core of LV finite growth models established by many authors [11][12][13][14][15].
One of the remarkable studies in this context was done by Göktepe et al. [11,12], which developed an efficient and a modular framework in order to perform finite element (FE) analyses of an athlete's heart, cardiac dilation and cardiac wall thickening.The myocardium was considered as a purely elastic isotropic Neo-Hookean material without the active response.The growth part of the deformation gradient was considered as an internal variable and prescribed according to the nature of the growth.For instance, they formulated myocardial growth in the athlete's heart as a pressure-driven, isotropic and volumetric growth process.On the other hand, the cardiac dilatation was characterized as strain-driven transversely isotropic.The total stretch along the myofibre direction was multiplicatively decomposed into an elastic and a growth part, which has a nonlinear strain-driven evolution and requires local Newton-Raphson iterations.Then, the growth part of the deformation tensor was characterized by a single scalar growth stretch along the fibre direction.The cardiac wall-thickening was, however, considered as stress-driven transversely isotropic growth.The growth part of the deformation gradient was determined through a single scalar growth stretch along the sheet direction.The framework was shown to be efficient and able to mimic the basic characteristics of the three growth types.
Kerckhoffs et al. [14] modelled both eccentric and concentric growth of the myocardium with a single strain-driven growth law, which was motivated by the statement that a strain-based stimulus correlates better with growth [16,17].It is assumed that both types of growth, i.e. eccentric and concentric, should be activated under volume and pressure overload and only the deformation state at end-systole and end-diastole should provide the necessary mechanical stimulus for growth, where the myocardium was formulated as a transversely isotropic elastic material and the electrophysiology was not considered [14].The solution procedure has basically three phases.They firstly determined the ventricle pressures via a non-linear lumped parameter system model of the systemic and pulmonary circulation [18].Then, the circulation stopped and the ventricles were allowed to grow.Once the growth was completed, the circulation was again started and implications of mitral regurgitation and aortic stenosis were studied.In other words, myocardial contraction and growth occur in three separate events at different times and they are all transiently intercorrelated.
Moreover, the eccentric growth was modelled in Lee et al. [15] by extending the framework of Göktepe et al. [12] to cardiac electromechanics in the sense of Sundnes et al. [19].The used growth stimulus was strain-driven and based on the average stretch that was calculated at the previous cycle.Since the growth tensor established using the previous cycle's stretch, it was treated as constant throughout the current cycle, which sequentially causes the derivative of the growth parameter with respect to the deformation field to be irrelevant.
More recently, Lamm et al. [20] developed a novel approach to describe stress-driven volumetric growth, which is capable of simulating both isotropic and anisotropic growth behaviour without using a predetermined growth tensor.Thanks to the not-predetermined growth tensor, heuristic assumptions are not required and the myocardium is not limited to grow along certain directions.Furthermore, the free energy function in the model consists of two components: the elastically stored energy in terms of the elastic right Cauchy Green tensor and growth-related energy in terms of the growth-related left Cauchy Green tensor.The growth component of Helmholz energy function is responsible for the description of cell expansion or contraction directly resulting in an accumulation of the growth-related energy due to the loading of the associated spring element.
Furthermore, many studies on myocardial growth have been conducted.We refer to Niestrawska et al. [21] for the overview of the current state-of-the-art of computational models of cardiac growth and remodelling.However, myocardial growth and ventricular remodelling have not been studied in an actively contracting patient-specific human ventricular model.Also, although the existence of viscosity in the myocardium has already been well accepted experimentally [22,23] and computationally [24,25], the viscosity of the myocardium has never been incorporated in finite myocaridal growth.In addition, in most of the previous studies, the myocardium was predetermined to grow in only one direction.However, in our model, the myocardium can experience the stretch-driven longitudinal (fibre) growth and stress-driven transverse (cross-fibre) growth, in order to imitate the serial and parallel addition of myocardial sarcomere.
The goal of the manuscript at hand is to analyze the effect of myocardial ischemia on ventricular remodelling and to evaluate the model's capacity to predict the pattern of LV remodelling under hemodynamic load with an integrated model of myocardial growth which embraces the active response and important mechanical features of the viscoelas-tic myocardium in continuously contracting patient-specific heart model.
For this purpose, the following examples are presented: • In the first example, we simulate acute regional myocardial ischemia in a patient-specific biventricular heart model and investigate the role of each growth type contributing to the ventricular remodelling in the ischemic region, where electrocardiogram (ECG) and the volumetime (v-t) curve are simultaneously computed.• In the second example, we compare the divergent patterns of LV remodelling when the heart is exposed to different pressure loading conditions and investigate its mechanism behind using LV heart model.
In finite element simulations for the myocardial growth, we use a fully implicit staggered solution procedure for the strong electromechanical coupling in a monodomain setting, where the myocardium is considered as an orthotropic electro-visco-elastic-growth material.Moreover, a standard Q1 element formulation is employed.The electrophysiological characteristics of the myocardium is modelled using the well-known phenomenological approach of Aliev and Panfilov [26].Furthermore, the ventricular blood pressure is coupled to the ventricle cavity volume and applied as a deformation dependent perpendicular load to the surface elements on the endocardium [27].The outline of this manuscript is as follows.Section 2 describes the kinematics, the balance equations and the constitutive equations for the electro-visco-elastic-growth behaviour in the cardiac tissue.Section 3 addresses the consistent linearization of the constitutive equations in a nonlinear finite element setting and the underlying local and global Newton iteration.In Sect.4, the constitutive equations along with evolution equations for the different types of growth and the type-selecting algorithm are introduced.In Sect.5, representative numerical examples that study the effect of myocardial growth on pathological ventricular remodelling are shown.In Appendix, we explain that the constitutive equations for the electrical and mechanical behaviour of the myocardium: cardiac electrophysiology, electrical flux, myocardial stress response and ventricular pressure.

Multiphysics problem: electro-visco-elastic-growth response
In this section, we introduce the fundamental variables and kinematics for the electro-visco-elastic-growth behaviour.Furthermore, the coupled field equations governing the myocardial growth, deformation and electrophysiology in the monodomain setting are explained.

Primary field variables and kinematics of the multi-field problem
The electromechanical state of the myocardium is described by the deformation map ϕ(X, t) and the transmembrane potential Φ(X, t) on the global scale.Moreover, the micromotion, n , λ a f }, due to the viscoactive deformation and the micro-diffusion, I g (X, t) = {λ g f , λ g s , λ g n }, due to the orthotropic myocardial growth, I e (X, t) = {c, r }, due to ion concentration alterations (e.g. [K] + , [Na] + , [Ca] 2+ ) across the cell membrane are treated as local fields.
The growing cardiac tissue is modelled in a continuum body B ∈ R 3 and B t ∈ R 3 at initial time t 0 ∈ R and current time t ∈ R, respectively.The macroscopic motion of the myocardium is described in the context of large strains by the non-linear deformation map that transforms the material points X ∈ B onto spatial points x = ϕ(X, t) at time t.The unit tangent of the Lagrangian point is mapped onto its counterpart in the Eulerian configuration by means of the deformation gradient F := ∇ X ϕ(X, t).The constraint J = det F > 0 guarantees non-penetrable deformations.Moreover, we utilize the symmetric right Cauchy-Green tensor and left Cauchy-Green tensor as a strain measure, respectively, where the covariant Eulerian metric tensor g ab = δ ab and contravariant Lagrangian metric tensor G AB = δ AB are used for index raising and lowering purposes in the Cartesian framework.The deformation gradient F is multiplicatively split into the mechanical part F ma and the growth part F g , where the mechanical part is further decomposed to elastic F e , viscous F v , and active F a contributions The viscous and the active as well as the growth part of the deformation gradient are prescribed along the three orthogonal directions (fibre f, sheet s and normal n), that are mutually orthogonal in the reference configuration.This simplification allows us to utilize the simplified and uncoupled one-dimensional evolution equations, which can be easily updated in a time-discrete setting.The active part of the deformation gradient is formulated along the fibre direction in terms of the active stretch λ a f , which is governed by the transmembrane potential Φ.Moreover, the viscous part of the deformation gradient with the orthotropic structure is defined as with viscous stretches λ v f , λ v s and λ v n related to three mutually orthogonal orientations f 0 , s 0 and n 0 , respectively.We use the rotation free viscous and active part of the deformation gradient that establish a fictitious intermediate state between the reference and current configuration in the sense of [28].Next, the growth part of the deformation gradient is formulated along the three directions where λ g f , λ g s and λ g n are the growth stretches along fiber, sheet and normal direction, respectively.The set of internal variables θ g = {λ g f , λ g s , λ g n } is introduced.The growth stretches are equal to one without growing features of the cardiac tissue.In case of the forward growth, the internal variables are larger than one, while in case of reverse growth the variables become less than one.Moreover, since the active, viscous and growth part of the deformation tensors are diagonal and rotation free as shown in Eqs. ( 3), ( 4) and ( 5), the order of F a , F v and F g in Eq. (2) becomes indifferent.Using the push-forward operation of the Lagrangian orientation vectors, the stretches in the three mutually orthogonal directions are obtained as Due to the orthogonal structure of the viscous, active and growth part of the deformation gradient, the stretches are further multiplicatively decomposed as Along fiber direction, the total stretch consists of elastic, viscous, active and growth parts.The total stretches along sheet and normal direction are decomposed into elastic, viscous and growth parts.These orthonormal stretches are shown in Fig. 1 along with the rheology for the electro-viscoelastic-growth response of the myocardium.Furthermore, the multiplicative form of the stretches along the orthogonal directions can be written in an additive format in the logarithmic strain space as follows

Balance equations
In what follows, /dt symbolize, respectively, the gradient and the divergence with respect to spatial coordinates x and the material time derivative of a quantity throughout the manuscript.Furthermore, the hat sign [ •] is employed to point out the variables that depend on the primary fields.The excitation-contraction coupling in the cardiac tissue is governed by two balance equations: the excitation equation in the monodomain setting and the balance of linear momentum, respectively, Herein, b depicts the volume-specific body forces in the reference configuration.Furthermore, the flux vector q, the ionic current Fφ and the Kirchhoff stress tensor τ have the following dependencies q = q(ϕ, Φ), F φ = Fφ (ϕ, Φ; I e ) and τ = τ (ϕ, Φ; I va , I e , I g ) .(10) Moreover, using the principles of thermodynamics, we can write the dissipation inequality as follows in terms of second Piola Kirchhoff S, Kirchhoff stress τ , the Helmholtz free energy ψ and the spatial rate of deformation tensor d, respectively.Additionally, the extra entropy term S 0 in the undeformed configuration is introduced in order to satisfy the second law of thermodynamics.

Constitutive equations of finite growth
In order to implement the growing viscoelastic myocardium, we adopt modified Hill model [29] the orthotropic viscoelastic response of the myocaridum that can be expressed in terms of the right Cauchy Green tensor of the intermediate configuration C ma = F − CF g−1 for the mechanical-active part, and the corresponding part of the deformation gradient F ma = FF g−1 .Following general theory of thermodynamics, the Helmholtz free energy ψ is introduced as and evaluate the dissipation inequality (11) σ σ σσ Fig. 1 a) Rheology for electro-visco-elastic-growth response of the myocardium in fibre direction, f. b) Rheology for visco-elastic-growth response of the myocardium in sheet and normal direction, s and n.
Both a and b are composed of the mechanical part and the growth part in serial.Furthermore, the mechanical part has two components in parallel.While the upper branch is the hyperelastic spring which represents the passive response of the myocardium in all the three orthogonal directions, the lower branch varies with the myocardial directions.In the case of the fibre direction f, it represents the electro-visco-elastic response, formulated through the elastic spring+dashpot+contractile element, but in the case of sheet and normal directions s and n, it represent elastic spring+dashpot Then, we can introduce the second Piola Kirchhoff stress tensor S as the conjugate to the total right Cauchy Green tensor C and the Kirchhoff stress τ as the conjugate to the covariant metric tensor g as follows where S ma is the mechanical-active second Piola Kirchhoff stress tensor, C ma the mechanical-active right Cauchy tensor and the mechanical-active Mandel stress tensor M ma = C ma S ma which is conjugate to the growth velocity gradient L g in the intermediate configuration.Furthermore, the mechanical-active Lagrangian moduli L ma and its Eulerian counterpart e ma are obtained by taking the second derivative of the free energy ψ in terms of C ma and g as follows Moreover, we introduce the evolution of the internal variables θ g , the diagonal components of F g , described based on deviations of a mechanically driven growth criteria φ g .
The set of growth stretches is determined based on the deformation gradient F at the end of the previous time step t n .We introduce the evolution equations for stretch-driven fibre growth and stress-driven sheet growth where φ g (F ma ) is the growth criterion for the case of stretchdriven longitudinal growth as a function of the mechanicalactive deformation F ma and φ g (M ma ) is the growth criterion for the case of stress-driven parallel growth as a function of the mechanical-active Mandel stress M ma = C ma S ma .
Remark: For the myocardial stress response, we adopt the modified Hill model [24] that describes the orthotropic electro-visco-elastic response, which is referred to as A.3.Moreover, the local dissipation inequality for the dissipative viscous damper is satisfied, which is presented in A.3.2.

Discrete equation of finite growth
In this section, we introduce the governing equations and their iterations for growth in the local and global level.We first present the local update of growth multipliers, where a finite difference scheme is used to discretize the evolution equation.Secondly, we introduce the nonlinear finite element method, which discretizes the equations for finite growth in space.

Local Newton iteration for finite growth
For the evolution of the growth stretches in Eq. ( 16), we use the limiting growth function k g (θ g ) proposed and validated in [11,30,31], and is written as where τ u g , τ d g , γ u g and γ d g are the adaptation speed parameters and (the degree of non-linearity of sarcomere deposition) for forward and reverse growth, respectively, and the subscript i = f, s, n, designates the orthogonal directions, fibre, sheet and normal in the myocardium.θ max and θ min are the sets of upper and lower limit for the growth stretches.To this end, we use the evolution of the growth stretches using the following finite difference approximation The evolution of the growth in Eq. ( 16) is approximated to formulate the discrete residual R through the utility of the finite difference approximation as follows where its sensitivities can be written as The set of growth stretches θ g is updated by using the following iteration

Global Newton iteration for finite growth
After the update of growth stretches through the local Newton iteration Eq. ( 20), the growth tensor F g and the mechanical-active deformation gradient tensor F ma in Eq.
(2) are obtained.Moreover, the second Piola Kirchhoff stress tensor S and the Kirchhoff stress τ are consecutively updated as presented in Eq. ( 14).Next, we introduce the fourth order tensor of Lagrangian moduli L as the derivative of the second Piola Kirchchoff stress tensor S in Eq. ( 14) 1 with respect to the total right Cauchy Green tensor C, which consists of four terms.The first term (2∂S/∂C) and second term (2∂S/∂F g ) are generic, while the third term (2∂F g /∂θ g ) and fourth term (2∂θ g /∂C) are derived differently and introduced in Sect.3, depending on the type of growth tensor.The first term is derived as and the second term is derived as where the different dyadic operators are defined as

Governing equations for myocardial growth
In the section, we introduce the governing equations for two different types of growth.Fundamentally, two types of mechanical stimuli exist that induce pathological LV remodelling: pressure overload and volume overload [32].The former leads to growth in myocardial thickness through additional sarcomeres in parallel (along sheet direction) while the latter cause growth in myocardial length via additional sarcomeres in series (along fibre direction).The simulation of these two different types of growth are stretch-driven fibre growth and stress-driven sheet growth in the myocardium, which are associated to the update of the growth stretches λ g f and λ g s , respectively.

Stretch-driven fibre growth
The stretch-driven fibre growth originates from the increased length of the sarcomere and the process can be expressed via a single growth stretch λ g f along fibre direction f 0 .The variables with subscript n are understood to be evaluated at time t n , and the variables without subscript are considered to be computed at time t n+1 .The total stretch along fibre direction is calculated as λ f = f 0 • F F • f 0 and is split as where λ ma f = λ e f λ v f λ a f and λ g f is the growth stretch along fibre direction.The growth component of the deformation gradient in Eq. ( 5) can be defined in terms of the growth stretch along fibre λ g f while other growth stretches are constant.Therefore, the derivative of the growth tensor with respect to the growth stretch along fibre λ g f at the time step the element undergoes the stretch-driven fibre growth is written as which is also the third term of the Lagrangian moduli in Eq. ( 21).The evolution equation for stretch-driven fibre growth is introduced in terms of the growth stretch along fibre λ where the evolution function for stretch-driven fibre growth k g and its derivative with respect to the growth stretch ∂k g /∂λ g f is written as where τ u f and τ d f are the time parameter for forward and reverse growth, γ u f and γ d f are the non-linearity parameter.Moreover λ max f and λ min f are the maximum and minimum growth stretch, respectively.The criterion for deciding whether to be forward or reverse growth along fibre direction is The residual explained in Eq. ( 18) for stretch-driven fibre growth is computed using an implicit Euler backward scheme ]/ t in order to secure unconditional stability as The local tangent K is defined in its specific expressions Then, the growth stretch along fibre direction is updated With the growth tensor at hand, F ma and S ma are obtained successively.Afterwards, the total stress tensor S can be calculated as in Eq. ( 14).With the tensors at hand, we can compute the fourth term in Eq. ( 21) as follows

Stress-driven sheet growth
Sarcomeres replicate along cross-fibre direction, leading to increased thickness of the wall in the heart, which is known as concentric remodelling [7].The governing equations for stress-driven myocardial growth along sheet direction is introduced in this section.Under stress-driven sheet growth, the growth stretch in sheet direction s 0 is only updated as λ g s ← λ g s,n and the other two growth stretches are not updated.The total stretch along sheet direction is calculated as λ s = s 0 • F F • s 0 and is decomposed as Its derivative with respect to the growth stretch along sheet λ g s is written as which is the third term of Eq. ( 21).The evolution equation for stress-driven sheet growth is introduced as follows where the growth stretch along sheet direction is determined using the following limiting function Unlike the growth along the fibre direction, where the criterion for growth is the total stretch along the fibre λ f , the reason for setting the term related to stress, mechanicalactive Mandel stress M ma as the criterion for growth along sheet direction λ s is to obtain more realistic simulations of the myocardium under the increased afterload.Afterload is defined as the ventricular wall stress, not stretch, that develops during systolic ejection.Stress-driven sheet growth is triggered when the trace of the Mandel stress of the mechanical-active part, tr(M ma ), exceeds the upward limit or falls below the downward limit where M u crit and M d crit are upper and lower limits of nongrowth stress range obtained during the healthy contraction.In case that tr(M) exceeds the upper threshold M u crit the forward growth along sheet occurs, while it falls behind M d crit reverse growth takes place in the sheet direction.The partial derivative φ g with respect to the growth stretch λ g s takes the following form where the sensitivities are computed as ∂C ∂λ ∂S ma ∂λ Moreover, we use an implicit Euler backward scheme with ]/ t so that the residual is computed in discrete forms as Using Eqs.(36) and (38), the derivative of the residual R is expressed as Then, the growth stretch along sheet direction is updated as With the updated growth stretch λ s , its derivative, which is also the fourth term in Eq. ( 21), can be obtained as where we calculate the term L 0 by pull back of the mechanical-active moduli L ma to the initial configuration

Type-selection algorithm
In this section, we introduce the algorithm for the selection of the growth-type in each element.The myocardium can experience stretch-driven fibre growth and stress-driven sheet growth, depending on its mechanical state.Moreover, forward or reverse growth can occur in both directions.We choose the gaps between the mechanical-active fibre stretch and its upper and lower limits, i.e., λ ma f −λ u crit or λ ma f −λ d crit , as the driving force for growth along myocardial fibre direction.We choose the gaps between the trace of the Mandel stress of the mechanical-active part and its upper or lower thresholds, i.e., tr(M ma ) − M u crit or tr(M ma ) − M d crit , as the drivers for growth along myocardial sheet direction.If the values of the mechanical-active fibre stretch λ ma f and the trace of the Mandel stress of the mechanical-active part tr(M ma ) neither exceed nor fall behind the upper and lower thresholds, in other words, the values are within the non-growth range, (λ d crit < λ ma f < λ u crit and M d crit < tr(M ma ) < M u crit ), then the both growth stretches are not updated.However, in case that the positive gaps occurs (φ g > 0), the forward myocardial growth is triggered so that the growth stretches along fibre and sheet increase, while the negative gaps (φ g < 0), cause the reverse myocardial growth so that the growth stretches along fibre and sheet become smaller than their previous values.Table 1 shows the algorithmic treatment of the selection of growth-type, stretch-driven fibre growth and stress-driven sheet growth as well as the computation of stress tensor and moduli.We refer to [29,33] for the details of model generation elling and hemodynamic disturbance induced remodelling, in order to show the validity of the introduced finite element framework.The human heart geometries, transmurally changing fibre direction, and boundary conditions are depicted in Fig. 2. The fibre angle in both heart models is linearly interpolated from 70 • on the endocardium to −70 • on the epicardium.The biventricular heart model is constrained by attaching linear springs to the nodes at the basal and epicardial surface with different stiffness values.For the LV heart model, the nodes at the basal surfaces are fixed and the nodes at the epicardial surfaces are constrained by linear springs with different stiffness values.In the first example, using the LV geometry shown in Fig. 2a, we simulate LV remodelling induced by two different hemodynamic disturbances; pressure overload and volume overload.In the second example, we simulate regional acute ischemia in the biventricular heart model shown in Fig. 2b, and investigate how the mechanical change due to ischemia eventually leads to ventricular remodelling, particularly focusing on which type of growth contributes to the remodelling.In both examples, the myocardial tissue is assumed to be in a resting state Φ 0 = −80 mV at time t = 0 ms.Unless differently specified, the required parameters for the examples are given in Table 2.Note that the ECGs in all simulations are recorded along lead II [34] by the integration of the electrical flux over the domain in the current configuration, i.e.B t q dv.Moreover, the scaling factor β t is formulated in terms of the activation time t act , which ensures that early depolarized elements have longer action potential duration (APD) than later excited elements.Moreover, the scaling factor β t is updated in every cycle and t act is the elapsed time measured in each element from the time when the atrioventricular node is stimulated to the time until the element reaches the transmembrane potential −40 mV.Before the final results are obtained, we perform many cycles until the electromechanical state of the material points become almost the same with the previous cycle so that we can obtain the ECG, the v-t curve, pressure-volume (PV) loops, which are almost the same with those obtained in the previous cycle.Moreover, we store the maximum and the minimum values of each elemental stretch and stress during the healthy saturated cycles, which will be used for the range in which the element does not grow.In other words, some disturbances leading the stretch and stress in elements beyond the normal range are necessary for the elements to initiate to grow.For the solution of the coupled balance equations, a staggered solution procedure is adopted as suggested by [33].
A constant time step t = 2 ms is used in all simulations.

Hemodynamic overload induced LV remodelling
In this section, we present a finite element analysis of an LV model, analyzing the influence of myocardial growth on the function, the remodelling and the mechanical behavior of LV in response to the hemodynamic overload.To do so, we design and execute two different simulations to consider two pathological cases, an LV with growing myocardium under   2. To imitate the volume overloaded LV, the pressure switch for preload p 1 is raised by 2 mmHg every cycle.Therefore, the preload in the LV, which is 4 mmHg at the beginning, reaches 14 mmHg in 5 cycles.In case of the pressure overload, the pressure switch for afterload p 2 is raised by 7 mmHg every cycle.Therefore, it increases gradually from 67 mmHg to 102 mmHg after 5 cycles.The snapshots of the simulations are shown in Figs. 3  and 4, which show the growth stretch along fibre direction λ g f and along sheet direction λ g s , respectively, obtained from LV with pressure overload and with volume overload.Prior to the application of the hemodynamic disturbance, the LV regularly beats with p 1 = 4 mmHg and p 2 = 67 mmHg, which does not show any growth along both directions.From the second cycle, we can observe the different pathological LV remodellings due to the disturbances in the LV inner pressure.The myocardium throughout the LV can grow along fibre direction or along sheet direction, as the type of growth in each element is selected at each time step, based on the mechanical environment.The typical features of concentric LV remodelling, such as wall thickening without the significant change of LV cavity volume [35], are manifested in the pressure overloaded LV, whose snapshots are in the left side of Figs.3a and 4a.Moreover, by comparing Figs.3a and 4a), it can be found that the extent of the growth along sheet direction λ g s is longer than that of the growth along fibre direction λ g f .Therefore, the growth type that contributes more to the concentric LV remodelling is the growth along sheet direction λ g s rather than that along fibre direction λ g f .Moreover, λ g s is more clearly visible in the mid LV region than in the basal and apical LV region, according to Fig. 3a and 4a.This is because the increase of the stress due to the LV pressure overload is revealed more strongly in the LV mid part, as the growth along sheet direction λ g s is triggered when the myocardium undergoes higher stress than the normal range.Furthermore, comparing the snapshots of the LV cross-sections, we notice that there is no substantial dissimilarities in the sphericity as cycles are repeated, despite the meaningful growth in the LV wall thickness.
The volume overloaded LV with increasing preload, on the other hand, develops a different pattern of remodelling.The wall thickness in the volume overloaded LV does not increase, but rather decreases, unlike the change of the wall thickness in the pressure overloaded LV.The wall thickness, measured at the longitudinal cross-section (z = −45 mm), is 10.01 mm in t = 800 ms prior to the volume overload of the LV and 8.95 mm in t = 4800 ms after 5 cycles under the volume overload.Moreover, we observe the progressively enlarged LV cavity with the cardiac dilation as cycles are repeated.These phenomena are accord with the characteristics of eccentric LV remodelling in the literature [14,36].Furthermore, comparing Fig. 3b and 4b, we found that the myocardial growth along fibre λ g f has a greater influence on the eccentric LV remodelling, rather than the growth along sheet direction λ g s has, which causes the initially elliptical shape of the LV gradually to be a circular shape.
The common characteristics of concentric and eccentric remodelling fundamentally originate from the myocardial fibre arrangement in each element over the LV heart model, as shown in Fig. 2, where we assign the fibre f with transmurally varying orientation from +70 • on the endocardium to −70 • on the epicardium as mentioned earlier.The vectors along myocardial sheet direction s are alinged in the transvere plane orthogonal to the fibre direction.To describe more clearly, we present Fig. 5 in which the vectors are projected on the representative planes.While the myocardial fibres are aligned along the LV, the vectors along the myocardial sheet directions are oriented in the radial direction.Therefore, the growth in the myocardial fibre stretch makes the LV shape more spherical, where the end-diastolic volume (EDV) increases from 198 ml to 293 ml after 5 cycles.On the contrary, the growth in the myocardial stretch along the sheet direction causes LV wall to be thicker without inducing the spherical shape and there is not any significant change in the inner LV volume, where the EDV decreases from 198 ml to 191 ml over 5 cycles.
Figure 6 presents the varied shape of PV curves accompanied by these two divergent patterns of remodelling.Figure 6a shows the PV loop obtained from pressure overloaded LV and Fig. 6b shows the PV loop from volume overloaded LV.In the pressure overloaded LV with the increasing afterload switch p 2 , the phase switching from isovolumetric contraction phase to ejection phase occurs at increasingly higher places in the PV loops.In addition, the highest pressure points, occurring during ejection phase, rises gradually, which means that the higher pressure required by the LV in order to pump blood throughout the whole body in each cycle.At the end, the myocardium undergoes higher stress than the stress occurred in the steady state with the normal afterload of 67 mmHg, which consequently causes stress-driven growth along sheet direction.The end-systolic volume (ESV), which is also considered to be the volume of residual blood remaining in the LV after ejection, gradually rises because it pumps less blood against stronger resistance due to afterload [37].Moreover, the volume of the inner LV cavity decreases as the sheet direction growth occurs which gradually reduces EDV.Due to the increasing preload value, the pressure applied to LV during filling becomes higher than the pressure during the filling without increased preload.Furthermore, at the end of the filling period, EDV increases because the LV is filled with higher blood pressure, and this rise in EDV indicates that the myocardium experiences a higher stretch than the stretch experienced in the normal state with p 1 = 4 mmHg.As a result, this increased stretch induces stretch-driven growth along fibre direction across the LV, which leads to the increased sphericity and the circular shape.Note that we do not conduct a comparative study between the simulations and the experimental results.We try to validate the applicability of our FE based framework through a personalized LV, which is differently remodelled in the end, depending on the type of hemodynamic overload.The personalized LV under pressure overload exhibits the increase in ventricular wall thickness with a little change in the LV chamber size caused by the myocardial growth along the sheet direction, which are accepted as the major features of concentric LV remodelling.In contrast, the LV under volume overload manifests the increase in the LV chamber size with a little change in wall thickness, which are known as the key characteristics of eccentric LV remodelling.These observations are well consistent to the literature [2,7,38,39].

Ischemia-induced remodelling in biventricular heart model
It is a known fact that the mechanical alterations due to myocardial ischemia can cause ventricular remodelling [5,40,41].In this example, we validate our FE framework for myocardial growth and ventricular remodelling due to regional acute ischemia in the biventricular heart model which is shown in Fig. 2b.After a healthy heart beat, acute regional ischemia is applied to the region of the heart model, which is shown in Fig. 7. To imitate the mechanical and electrical changes induced by the ischemia, the following parameters related to the passive tissue stiffness (κ, a, a f , a s , a f s ) are decreased to 25% of the normal values [42,43], the contractility parameter (ξ ) are decreased to 10% of the normal value, the electrical conduction parameters d iso and d ani are decreased to 25%, and the parameter for Fig. 3 The snapshots of the distribution of growth stretch in the fibre direction, taken at end-diastole.a) illustrates the snapshots of the LV under pressure overload where the afterload is increased by 7 mmHg every cycle, finally reaching p 2 = 102 mmHg.b) depicts the snapshots of the LV under volume overload where the preload is increased by 2 mmHg every cycle, until it reaches to p 1 = 14 mmHg.The volume overloaded LV heart model becomes more spherical, while the pressure overloaded one experiences wall thickening.It is in the volume overloaded LV heart model that the growth stretch in the fibre direction is more prominent, not in the pressure overloaded one the action potential duration γ is increased to 125% of the normal value [44].The parameters are adjusted to new values linearly in four cycles (in 3200 ms).Moreover, the evolutions of the internal variables λ g f and λ g s in the ischemic regions are shown in Fig. 8, and the snapshots obtained during simulation are presented in Fig. 9 with ECG and v-t curve recorded throughout the cardiac cycles.The snapshots presented in Fig. 9 are taken at EDV state every 2 cycles, immediately before the onset of contraction.Figures 9a and 9b   The decreasing contractility and stiffness due to ischemia lead to an unbalance between the contractile forces of the ischemic segments and the normal segments.As the cycles repeat, the myocardium in the ischemic regions is more stretched, which stimulates stretch-driven fibre growth in the myocardium and develops the larger myofibre stretches in the region.Subsequently, we could observe that the local ischemia region becomes more spherical and thinned, losing the ellipsoid shape of the ventricle with the increased growth stretch along fibre direction λ Fig. 6 PV loops recorded during the simulations of the pressure overloaded and volume overloaded LV heart model: a)The simulation of the pressure overloaded LV with increasing afterload.As the switch value p 2 increases, the internal pressure that the LV must overcome in order to switch the phase from isovolumetric contraction to ejection gradually increases.Therefore, the maximum LV pressure also increases each cycle.The increase of the EDV is not prominent but, the increase of the ESV stands out.b)The simulation of the volume overloaded LV with increasing preload.As the switch value p 1 increases, the internal LV pressure during filling phase increases as well.Consequently, the LV heart model is progressively dilated with the increasing EDV and ESV each cycle Fig.7 Illustration of the ischemic zone from different views, in which the ischemic regions are shown in darker grey.In order to mimic mechanical and electrical changes induced by ischemia in the regions, the parameter for contractility (ξ ), for stiffness (κ, a, a f , a s , a fs ), for the conductance velocity (d iso and d ort ), and for AP duration (γ ) are modified with the characteristic of the remodelling due to the myocardial ischemia in the literature [41] The amount of the increasing stretch-driven fibre growth over the region is not identical despite the fact that the reduction in stiffness is identical throughout the ischemic region.The higher extent of the growth is observed in the center of the ischemic region than in the border as shown in Fig. 9. Furthermore, the myocardium in the ischemic region undergoes the reverse growth along sheet direction due to less contractility of the myocardium induced by the ischemia.The passive stress in the less contracting myocardium gradually decreases and tr(M), which is the cause of the stress-driven sheet growth in the myocardium, falls below the lower limit obtained in the healthy state, which generates the reverse growth along sheet direction with λ g s < 1 [-].In other words, the reverse myocardial growth in the sheet direction λ g s occurred by a decrease in stress measured by tr(M) due to the decreased contractility, which contributes to the wall thinning in the ischemic region.The fact, that the wall thinning is accompanied by the reduced wall stress, is a conflicting b)The snapshots taken at 7200 ms show the contour of the growth along fibre and sheet direction, and the locations of the nodes under investigation.Node A is closely located at the border of the ischemic region, and Node B and C are near the center of the region.We observe that the myocardial growth tends to be large nearby the center of the ischemic region, and decreases toward the periphery.However, there is a region where growth along sheet direction λ g s does not occur even in the center of the ischemic region result to the literature, where the increase of wall stress and wall thinning are observed together in the infarcted region in the ventricle [41].
Figure 8 depicts nodal values of λ g f and λ g s which are observed at the nodes A, B and C in the ischemic area.While node A is located closer to the border of the ischemic region, B and C are in the middle of the region.As mentioned earlier, the extent of myocardial growth is not identical over the region.The extent of growth is relatively lower at the border region, and the closer the node locates to the center of the region, the larger growth is observed.For example, we observe that the growth along fibre direction is the largest at node C, followed by B, and then A. A similar trend is observed in the growth in the sheet direction s, but there are exceptions.For example, the extent of growth along the sheet direction λ g s at node B, located in the middle of the ischemic region, is lower than the growth at node A, which is nearby the border region.In the outside of the ischemic region, almost no growth is present and even if it occurs, the degree of growth is very small.Moreover, it is worth to note that the both increasing and decreasing stretches are forming a step-like progression.If the step-jump in the growth along sheet direction λ g s occurs, then it is at the systolic phase.On the other hand, the rise of the growth along fibre λ g f is mostly generated at each diastolic phase.For example, the first increase of λ g s is developed around 870 ms because this is the first systole after the onset of regional ischemia, and the first step-jump in λ g f is produced at about 1550 ms because this is the first diastole after the initiation of ischemia.As both stretches for myocardial growth (λ g f and λ g s ) in elements accumulate, the regional distension and eccentric remodelling in the ischemic area are induced, which in turn leads the LV to be more spherical, as shown in Fig. 9.The right ventricle, except the effect from the ischemic region in septum, nearly remains unaffected from remodelling, as there is no ischemic region.Furthermore, Fig. 9c illustrates the corresponding ECG and v-t curve.The slowed conduction speed and the reduced AP duration in the ischemic region have influence on the electric domain.These effects are not significant enough to cause heart failure, but we can observe that the heart model is electrically disturbed, and it manifests the typical characteristics of ischemia in ECG, such as more peaked QRS complexes and the hightened and widened T wave [44].From the v-t curve, it is observed that the LV cavity volume steadily increases as the distention of the ischemic region progresses.EDV and ESV are also continuously rising, which consequently hinders the cardiac function in terms of EF.

Conclusion
In this work, we showed the applicability of a novel threedimensional material model that describes an electro-viscoelastic-growth formulation on the myocardium in a finite element based numerical framework by using a personalized LV and biventricle model.A fully implicit staggered solution procedure is utilized in the FE framework for the strong electromechanical coupling.Unconditional stability of time integration is ensured as we solved all time integrations through an implicit backward Euler integration scheme.
Unlike to other growth models in the literature where the growing direction of the myocardium is predetermined in one of the three orthogonal directions, the growth direction for each element is selected at each time step based on the mechanical conditions.Furthemore, the myocardium can be returned to the initial configuration or shrinked to a stretch value less than 1 [-].With this model, we investigated how the Fig. 9 Snapshots, ECG and v-t curve that demonstrate ventricular remodelling caused by acute regional ischemia.The contour of growth stretch along fibre λ g f is shown in a), the contour of growth stretch along sheet λ g s is shown in b), and the graphs in c) illustrate ECG where the normalized voltage (upper curve) is plotted on y-axis, and LV v-t diagram where the normalized volume (lower curve) is plotted.After the regular beat in the first cycle without ischemic region, the steady contracting biventricular heart model is disturbed by acute regional ischemia over the four-cycle period.As mechanical stiffness and the contractility decrease, the ischemic region undergoes higher stretches, which consequently causes stretch-driven growth along fibre direction as well as the decreased contractility produces less stress, which generates stressdriven reverse growth.c) shows the ECG that reveals the effect of the slowed conduction velocity and the reduced APD.ECG progressively exhibits taller QRS complexes and taller and more widened T wave.The increased LV cavity volume and the deteriorated cardiac function at each cycle are observed in the v-t curve, where the EF values are shown shape of the entire heart geometry changes due to the accumulation of myocardial growth at the element level, which agrees well with literature [2,7,38,39].
In the first example, pathological ventricular remodelling due to hemodynamic disturbances, pressure overload and volume overload, were simulated using the LV heart model.The hemodynamic alteration at organ level causes the myocardium to grow in different directions at the element level, which leads the LV heart model to the different shapes of remodelling.The pressure overloaded LV heart is simulated by increasing the afterload every cycle, which generates the higher peak pressure during the ejection period.Therefore, the elements throughout the LV experience the stress driven growth in the sheet direction more predominantly, rather than stretch-driven growth along fibre direction.As the sheet-direction growth in individual elements is aggregated macroscopically in the whole heart level, we observe finally LV wall thickening, which is one of the features of LV concentric remodelling.Meanwhile, the volume overloaded LV heart is simulated by increasing preload at every cycle.Exposed to a higher preload than the regular one, the LV is filled with a larger amount of blood volume, which increases the EDV, and consequently the myocardium is subjected to higher stretch than that of the normal state.This primarily produces stretch-driven growth in individual elements rather than stress-driven growth.Subsequently, this leads the LV model to undergo an eccentric remodelling, revealing its typical characteristics such as LV dilation with wall thinning as the cardiac cycles are repeated.
Secondly, we simulated acute regional ischemia in the biventricular heart model, in which the ischemic region experiences mechanical and electrical alterations in myocardial tissue.These disturbances in the ischemic tissue induce the myocardial fibre growth and reverse sheet growth, subsequently the region becomes progressively stretched and shows ventricular wall distention.At the element level, we observe the increased total fibre stretch due to stretch-driven fibre growth, which leads the ischemic region to be dilated, while the decreased contractility causing the myocardial stress to be lowered in the region contributes to the wall thinning through the reverse growth along the myocardial sheet direction.Conclusively, we observe that the LV shape becomes more spherical due to the myocardial growth and ventricular remodelling.Furthermore, the slowed conduction velocity and the shortened APD produce the gradually widened and hightened T wave in ECG as cardiac cycles are repeated.
However, there is some room for improvement in the framework.The influence of viscosity versus elasticity in the myocardial growth and remodelling was not studied, because we did not compare our simulations to the elastic cases.Nevertheless, the reason why we use the viscoelastic material model, rather than a purely elastic model, is that viscoelasticity is an indispensable aspect as mentioned in the literature, especially when studying the mechanical behaviour of the myocardium [22][23][24][25].Therefore, we plan to conduct a detailed comparative study on the effect of different material models on the myocardial growth.Moreover, in the second simulation in Sect.5.2, we did not consider the ventricular inner pressure.The mechanical deformation could be more significant if there is a reaction force of the pressure, which consequently affects the remodelling pattern.Furthermore, since mechano-electric feedback (MEF) was not considered in the simulations, the disturbance in cardiac electrophysiology was caused only by the alterations in APD and conduction velocity due to ischemia.As a result, the effect of electrical disturbance to the LV function was not significant.However, we expect that the electrical disturbances can be additionally caused by myocardial growths through MEF.More specifically, we hypothesize that the myocardial growth in the fibre or sheet direction does not only induces the mechanical changes, but also generates electricity through different types of stretch activated channels [45], which will eventually affect cardiac electrophysiology.We plan to simulate ischemia-induced premature ventricular contraction with growing myocardium under MEF in future studies.
where a, b, a f , b f , a s , b s , a fs and b fs are constants describing the deformation state of the isotropic and orthotropic characteristics of the myocardium in terms of the invariants Macaulay brackets x = (x +|x|)/2, which exclude the negative values of I 4i − 1, are adopted because the myocardium does not bear any compressive load [54].Moreover, we treat the passive myocardium as a quasi-incompressible material, which requires a generic free energy function for the volumetric part, which is employed as Then, the passive part of the Kirchhoff stress tensor is written as (A12) The passive part of the tangent moduli L p is additively decomposed into two contributions, L vol and L p , which are written as where the coefficients are (A14) The fourth order tensor I C ma −1 in equation (A13) is (A15) The fourth order structural tensors, F and S, are defined as F = f 0 ⊗f 0 ⊗f 0 ⊗f 0 and S = s 0 ⊗s 0 ⊗s 0 ⊗s 0 .(A16)

A.3.2 Visco-active stresses
In order to incorporate the visco-active features of the myocardium, a simple quadratic free energy function in the logarithm strain space is considered through (A18) Remark: The viscous part of free energy function in terms of logarithmic strain is not polyconvex.However, it can be shown that the rank-one convexity is satisfied for the logarithmic contribution of the free energy function, see Cansız et al. [24].

Evolution equations for internal variables
Based on the Biot equation [55], the following evolution equations for the viscous stretches are obtained as in terms of material parameters m i , η i and βi .Note that the activation stress βi is considered as unity to ensure the consistent units.Furthermore, the active stretch representing the myocardial contraction is modelled as the following exponential type evolution equation where k 1 , ζ and p are the material parameters.

Local dissipation inequality
In the sense of Biot [55], the non-equilibrium response of the material can be worked out from the relatation between a storage and dissipation potential where the non-equilibrium free energy function ψ v characterizes the energy stored in the spring of the Maxwell branch and the dissipated energy in the dashpot is described by the dissipation potential D. The viscous part of the free energy function is a simple quadratic form in the logarithmic strain space as presented in Eq. (A17).Furthermore, the relation between the viscous free energy function and the logarithmic stresses and thermodynamical forces conjugate to the internal variables are expressed as and For the dissipation potential, a power-type generic expression is utilized governed by the material parameters m i , η i and β i .Herein, m i describes the transition from purely linear viscoelastic to nonlinear viscoplastic flow, η i is the viscosity parameter and β i designates the activation stress and will be set as unity to ensure the consistent units.Moreover, an evolution equation for the internal variables can be attained by making use of Eqs.(A22-A24) and some manipulations, Eq.(A19) of the manuscript is obtained.Thanks to the definition of β i and the rate of the internal variable v i , the local dissipation inequality naturally is satisfied, i.e. β i v i ≥ 0.

A.4 Ventricular blood pressure
A diagram for PV relationship is described in Fig. 10.We assume the ventricular pressure is quasi-linearly dependent on the volume in isovolumetric contraction, isovolumetric relaxation, and filling.Therefore, a linear evolution law is adopted as ṗv = κ( θ − θ) with θ := Ṽ V 0 and θ := Therein, κ is the material parameter having different values for each phase (isovolumetric contraction, isovolumetric relaxation, and filling) in order to enforce the condition V ≈ 0 during isovolumetric phases and ṗv ≈ 0 during filling.Furthermore, θ , θ , Ṽ and V 0 represent the dilatation at the end of the preceding phase, the current dilatation, the volume at where the material parameters C ap , R c and R p determine the resistance and compliance properties of the blood flow and q represents the outflow from the ventricle.Figure 10 shows also the switch conditions from one phase to the next one.
The phase movement from filling to isovolumetric contraction takes place if q > 0, from isovolumetric contraction to ejection if p v > p 2 , from ejection to isovolumetric relaxation if q < 0 and from isovolumetric relaxation to filling if p v < p 1 .
Further information about the computation of the stress tensor and the moduli regarding the myocardial stress response of modified Hill model is explained in A.3.In addition, we present the Kirchhoff stress tensor τ which is obtained through the push forward of the second Piola Kirchhoff stress tensor S and the Eulerian moduli e as follows τ = FSF and e = [F⊗F] : L : [F ⊗F ] .(24)

14 goto 31 !
Fig. 2 a) Virtualized LV geometry discretized by 18809 four-node tetrahedral Q1 elements and 4116 nodes and b) virtualized biventricle geometry discretized by 25366 four-node tetrahedral Q1 elements and 6270 nodes.From left to right: the discretized geometry along with the dimensions (the lengths are in millimetres.),two different views of the fibre orientation and its angle θ distribution and the springs with

Fig. 4
Fig. 4 The snapshots of the distribution of growth stretch in the sheet direction, taken at end-diastole.a) illustrates the snapshots of the LV under pressure overload where the afterload is increased by 7 mmHg every cycle, finally reaching p 2 = 102 mmHg.b) depicts the snapshots of the LV under volume overload where the preload is increased by 2

Fig. 5 a
Fig. 5 a) and b) present the alignment of the myocardial fibre f and sheet s vectors in the representative cross sections, respectively.c) shows the locations of the cross sections in the LV heart model

Fig. 8 a
Fig. 8 a)The evolutions of nodal growth λ g f and λ g s in different regions.b)The snapshots taken at 7200 ms show the contour of the growth along fibre and sheet direction, and the locations of the nodes under investigation.Node A is closely located at the border of the ischemic region, and Node B and C are near the center of the region.We observe that the ψva (C ma , s; F e ) = the shear moduli μ i and elastic logarithmic strains e i = i − v i − a i in the mutually orthogonal directions of the myocardium, respectively.Then, the visco-active part of the Kirchhoff stress tensor is described asŜva = 2∂ I 4f ψf va f 0 ⊗ f 0 + 2∂ I 4s ψs va s 0 ⊗ s 0 + 2∂ I 4n ψn va n 0 ⊗ n 0 .

Table 1
Algorithm for the selection of the growth and the update of the stress tensor and

Table 2
Parameters used in the simulations Mechanics