Computational modeling of coupled cardiac electromechanics incorporating cardiac dysfunctions

European Journal of Mechanics A/Solids xxx (2014) 1e14

Contents lists available at ScienceDirect

European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol

Computational modeling of coupled cardiac electromechanics incorporating cardiac dysfunctions
lu, H. Onur Solmaz, Serdar Göktepe* Ezgi Berberog Department of Civil Engineering, Middle East Technical University, TR-06800 Ankara, Turkey

a r t i c l e i n f o

a b s t r a c t

Article history: Received 20 November 2013 Accepted 28 February 2014 Available online xxx

Computational models have huge potential to improve our understanding of the coupled biological, electrical, and mechanical underpinning mechanisms of cardiac function and diseases. This contribution is concerned with the computational modeling of different cardiac dysfunctions related to the excitation econtraction coupling in the heart. To this end, the coupled problem of cardiac electromechanics is formulated through the conservation of linear momentum equation and the excitation equation formulated in the Eulerian setting and solved monolithically through an entirely finite element-based implicit algorithm. To model the electromechanical coupling, we use the recently proposed, novel generalized Hill model that is based on the multiplicative decomposition of the deformation gradient into the active and passive parts and on the additive split of the free energy function. This framework enables us to combine the advantageous features of the active-stress and the active-strain models suggested in literature. The proposed coupled approach is further supplemented by the Windkesel-based model to account for the pressure evolution within the ventricular chambers during the cardiac cycle. This allows us to generate the pressureevolume curves as a diagnostic tool to detect possible cardiac dysfunctions and to assess the efficiency of the heart function. The proposed model is employed to investigate different pathological cases that include infarction, eccentric hypertrophy, and concentric hypertrophy. The effects of these distinct cardiac dysfunctions on the pressureevolume curves and on the overall excitationecontraction of the heart are computationally examined and compared to the clinical observations reported in literature.
2014 Elsevier Masson SAS. All rights reserved.

Keywords: Coupled cardiac electromechanics Cardiac diseases Pressureevolume curves

1. Introduction Heart disease is not only the most common cause of mortality in industrialized countries but also cause a huge amount of financial cost worldwide. More importantly, the rate of deaths due to cardiac disease is expected to rise in the near future (Go et al., 2013). In pathological cases, the cardiovascular system may become inefficient and dysfunctional in maintaining the vital functions of the body. The side effects of cardiovascular malfunctions may result in the disabilities of other organs as well. In case of a heart failure, the heart cannot supply the peripheral tissue and organs with a sufficient amount of blood flow resulting in a lack of oxygen delivery (Klabunde, 2011). This type of disorders may arise from chronic valve disease, chronic arterial disease, and the failure of the myocardium itself, (Formaggia et al., 2010). In order to reduce the

* Corresponding author. E-mail addresses: (S. Göktepe).

[email protected],

[email protected]

rate of death and to improve the quality of patients’ life, diagnostic techniques and treatment procedures should be easily accessible, time-saving, and as inexpensive as possible. Therefore, the understanding of underlying mechanisms controlling the heart function and diseases is of key importance for all steps of diagnostic and therapeutic techniques. With this overarching goal in mind, concurrent to the purely medical and biological research, the interdisciplinary studies involving cardiologists, cardiac surgeons, biologists, material scientists, and engineers have been intensively undertaken over the last decades. As it is the case for man-made materials, the predictive constitutive models furnished with sound and robust numerical tools for computational cardiology play a central role in guiding, improving, and optimizing the mainly trial-and-error based diagnostic and surgical techniques (Kerckhoffs et al., 2009; Constantino et al., 2012). Today, the electrocardiograms (EKG) and the pressureevolume (PV) curves are among the most broadly used non-invasive methods to diagnose heart diseases. The EKGs are especially used for the diagnosis of pathological cases related to the cardiac electrical conduction system that cover rhythm disturbances, blockages

http://dx.doi.org/10.1016/j.euromechsol.2014.02.021 0997-7538/
2014 Elsevier Masson SAS. All rights reserved.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

2

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

in the electrical conduction system, myocardial ischemia, and infarction. The PV curves describing the change of the ventricular pressure with respect to the ventricular volume in each cardiac cycle serve as a diagnostic tool to assess the efficiency and the work expended by the cardiovascular system in the case of blood flow abnormalities (Klabunde, 2011). In addition to their use as diagnostic tools in clinical studies, the EKGs and the PV curves can also be generated computationally. Therefore, they can be used as a metric to assess the predictive modeling capability of the computer models and serve as a bridge between clinical findings and computational simulations. In this work, we investigate the effects of different cardiac dysfunctions on the efficiency of the heart by using the PV curves. Hence, it is worth giving a brief introduction about the main characteristics of the PV curves. A schematic representation of the pressureevolume relationship of the left ventricle is depicted in Fig. 1. The cardiac cycle is initiated by the electrical impulses generated at the sinoatrial (SA) node. These electrical waves first spread throughout the atria and enter the ventricles through the atrioventricular (AV) node. These signals are then transmitted to the ventricular cardiac muscle cells (myocytes) by the fast conducting His bundles and the Purkinje fibers in sequence and result in the contraction of the ventricles upon depolarization of the myocytes. The myocytes then repolarize back to their resting state through the local transmembrane ion dynamics. This leads the ventricles to resume their relaxed state (Katz, 2011). The cardiac cycle can be divided into four main phases: isovolumetric contraction, ejection, isovolumetric relaxation, and filling (Fig. 1). The sequence of these four phases follows the pressure loop counter-clockwise. In Fig. 1, the cardiac cycle starts at Point A that corresponds to the end of filling phase (end-diastole). The first phase, between Points A and B, refers to isovolumetric contraction where the left-ventricular pressure (LVP) increases abruptly without a significant change in the left ventricular volume (LVV). In this phase, the left ventricle starts to contract while the mitral valve is closed. Thus, no blood outflow from the left ventricle takes place during isovolumetric contraction. As the ventricular

Fig. 1. Schematic illustration of the left ventricular pressure (LVP)-left ventricular volume (LVV) curve (solid line) demonstrating the four main phases of the cardiac cycle: isovolumetric contraction, ejection, isovolumetric relaxation, and filling. The end-systolic pressureevolume relationship (ESPVR) and the end-diastolic pressuree volume relationship (EDPVR) are depicted by dashed curves. The difference between the end-diastolic volume (EDV) and the end-systolic volume (ESV) is called the stroke volume (SV).

pressure exceeds the aortic pressure, the ejection phase, between Points B and C, begins at Point B. The ejection phase is accompanied by continuing contraction of ventricles and a decrease in the LVV due to the blood ejected through the aortic valve. The systolic phase starts at Point A and ends at Point C where the diastolic phase starts with the isovolumetric relaxation as the LVP drops below the aortic pressure. During the isovolumetric relaxation, the LVP continues to decrease and the LVV remains constant. As the LVP falls below the left atrial pressure at Point D, the mitral valve opens and the filling phase starts. Due to the blood flow from the atria, the LVV increases until Point A. The four phases are repeated for each cardiac cycle and the transition between these phases is controlled by the pressure gradients between the left ventricle, aorta, and left atrium. The right ventricular PV curves possess a similar shape and phases. Owing to the structural differences between the right and left ventricles, however, the magnitude of the right ventricular pressure is about one-fifth of the LVP. In the assessment of the ventricular pumping efficiency, the stroke volume (SV), the difference between the end-diastolic volume (EDV) and the end-systolic volume (ESV), is used as a basic indicative measure. In addition, the end-diastolic and the end-systolic pressureevolume relationships, denoted as EDPVR and ESPVR in Fig. 1, are associated with the stiffness of the overall heart in the most relaxed and contracted states, respectively (Burkhoff et al., 2005). Computational models serve as an important tool providing enormous potential to improve our understanding of the coupled bio-electro-mechanical underpinning mechanisms of cardiac function and cardiac diseases. The modeling of cardiac tissue, however, involves several challenges that are far from being resolved completely. The heart function is governed by concurrent electrical and mechanical activities occurring at different scales. Furthermore, these mechanisms mutually affect one another and thus require sound constitutive models and robust numerical techniques accounting for the strongly coupled electromechanical response. Moreover, cardiac tissue has an inherently anisotropic micro-structure (Fig. 5), vital for making the heart a threedimensional biological pump. Cardiac tissue is a nonlinear material that undergoes large deformations. Hence, the developed heart models must take the geometrical and material nonlinearities into account. Besides, the inner surface of the ventricles is subject to transient pressure due to the force exerted by blood during the course of every heart beat. The constitutive approaches to the electro-mechanically coupled response of the heart proposed in literature can be classified into two main groups: the active-stress models and the active-strain models. In the active-stress-based models, the total stress tensor is additively decomposed into the passive and active parts. While the passive part of the stress tensor is modeled as a function of deformation only, a transmembrane potentialdependent evolution equation is used to model the active part of the stress tensor, see Nash and Panfilov (2004), Panfilov et al. (2005), Nickerson et al. (2006), Keldermann et al. (2007), Niederer and Smith (2008), Göktepe and Kuhl (2010), Lafortune et al. (2012), Eriksson et al. (2013a,b). The active-strain-based approaches go back to the works by Nardinocchi and Teresi (2007) and Cherubini et al. (2008) where the deformation gradient is multiplicatively split into the active and elastic parts. The stress tensor is then obtained from an energy storage function, which depends only upon the elastic part of the deformation gradient. This approach has been extended by Ambrosi et al. (2011), Nobile et al. (2012), Rossi et al. (2012) toward anisotropic constitutive models for cardiac tissue and by Stålhand et al. (2011) to model smooth muscle contraction. These two approaches have different advantages and disadvantages one over another from different viewpoints. The additively decomposed stress response in the


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

active-stress approaches can be advantageous when the associated material parameters are identified by using different sets of experimental data acquired in the absence of electrical stimulation (passive) and during the active contraction, where the active and passive mechanisms co-exist. From the experimental point of view, however, the direct measurement of strains is easier than the measurement of the local stress fields, which requires indirect techniques. In addition, the tools of mathematical analysis is more readily available to judge the stability of a material model when the stress tensor is obtained from an energy storage function as it is the case in the active-strain models. One of the main motivations for the development of computational heart models is to develop new diagnostic tools to advance the existing therapeutic techniques (Kerckhoffs et al., 2006). Although the computational tools have been widely utilized to investigate electrophysiological disorders such as arrhythmias, bundle blocks, and fibrillation (Dal et al., 2012; Kerckhoffs et al., 2012; Rantner et al., 2013; Roberts et al., 2012), the computational models accounting for electromechanical cardiac dysfunctions are rather limited. Among these the reader is referred to the works by Kerckhoffs et al. (2007) and Kerckhoffs et al. (2012) where the PV curves have been utilized to investigate the ischemia, pulmonary artery constriction, and the blockage of conduction in the left bundle. In the present work, we develop computational heart models incorporating infarction, eccentric hypertrophy, and concentric hypertrophy. For this purpose, the coupled problem of cardiac electromechanics is formulated through the balance of linear momentum and the excitation equation formulated in the Eulerian setting. To model the electromechanical coupling, we use the recently proposed, novel generalized Hill model (Göktepe et al., 2013a,b) where the classical Hill model (Hill, 1938, 1970) has been extended to the three-dimensional finite deformation setting. Following the active-strain models, we decompose the total deformation gradient into the active and passive parts. The active part of the deformation gradient is considered to be dependent upon the intracellular calcium concentration whose evolution is driven by the action potential. Furthermore, the inherently anisotropic, active micro-structure of cardiac tissue is accounted for in the tensorial representation of the active part of the deformation gradient. As opposed to the active-strain models where merely the active-passive split of the deformation gradient has been utilized, we further additively decompose the free energy function into passive and active parts. This decomposition allows us to recover the additively split stress structure of the active-stress models. Therefore, the advantageous features of the approaches that employ either additive stress decomposition or multiplicative split of the deformation gradient to account for excitation-induced contraction are incorporated within a unified constitutive framework. The two-parameter Alieve Panfilov model (Aliev and Panfilov, 1996) is utilized to simulate the electrical excitation (Göktepe and Kuhl, 2009). The fully coupled cardiac electromechanics problem is solved through a unified implicit finite element approach. The algorithmic formulation of the geometrically non-linear, strongly coupled problem is carried out in the Eulerian setting resulting in unconditional stability and quadratic convergence (Göktepe and Kuhl, 2010). The proposed coupled approach is further supplemented by the Windkesel-based model to account for the transient pressure evolution within the ventricular chambers during the cardiac cycle. This enables us to generate the PV curves (Fig. 1). The proposed model is employed to investigate different pathological cases that include infarction, eccentric hypertrophy, and concentric hypertrophy. The distinct effects of these cardiac dysfunctions on the PV curves are examined by the model

3

developed. The predictive capacity of our model is demonstrated by comparing the computationally obtained results with the clinical findings. The manuscript is organized as follows. Section 2 addresses the three-dimensional continuous formulation of the coupled electromechanics of the heart. In Section 3, the specific forms of the constitutive equations are particularized. In Section 4, the fully coupled, three-dimensional finite element analysis of a generic heart model is conducted to investigate the effect of infarction, eccentric hypertrophy, and concentric hypertrophy on the PV curves. Finally, Section 5 concludes with some critical remarks. 2. Coupled cardiac electromechanics This section is concerned with the general three-dimensional formulation of the coupled electromechanics of the heart. A coupled initial boundary-value problem of cardiac electromechanics within the Eulerian setting is formulated in terms of two primary field variables: the placement 4(X, t) and the transmembrane potential F(X, t). The placement is the non-linear deformation map, depicted in Fig. 2, and the transmembrane potential is defined as a potential difference between the intracellular medium and the extracellular medium within the context of monodomain formulations of cardiac electrophysiology (Keener and Sneyd, 1998). The temporal and spatial evolution of the primary fields is governed by the balance of linear momentum and the reactionediffusion-type equation of excitation, which are introduced in Section 2.2.

2.1. Kinematics: activeepassive decomposition The reference configuration B3R3 and the current configuration S3R3 at time t˛Rþ of an excitable and deformable solid body are depicted in Fig. 2. The non-linear deformation x ¼ 4t ðXÞ : B/S maps the material points X˛B onto their spatial positions x˛S at time t. The deformation gradient, defined as F : ¼ VX 4t ðXÞ : TX B/Tx S, is the tangent map between the tangent spaces of the respective configurations. Here, the gradient operator Vx[] denotes the spatial derivative with respect to the reference coordinates X and the Jacobian J: ¼ det F > 0 is the volume map. Motivated by the kinematics of finite plasticity (Kröner, 1960; Lee, 1969) and the work by Cherubini et al. (2008), the deformation gradient is multiplicatively decomposed into the passive part Fe and the active part Fa

Fig. 2. Motion of an excitable and deformable solid body in the Euclidean space R3 through the non-linear deformation map 4t(X) at time t. The deformation gradient F ¼ VX4t(X) describes the tangent map between the respective tangent spaces. The deformation gradient F is multiplicatively decomposed into the passive part Fe and the active part Fa with Ba denoting the fictitious, incompatible intermediate vector space.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

4

The second partial differential equation of the coupled problem, the excitation equation,

i h f b  bI ¼ 0 _F  J div J 1 q

in

B

(6)

describes the spatio-temporal evolution of the action potential field b  and the non-linear F(X, t) in terms of the diffusion term div½J 1 q f _ : ¼ D½=Dt stands for the macurrent source bI . The notation ½ terial time derivative. Within the framework of FitzHughe Nagumo-type models of electrophysiology (Fitzhugh, 1961; Fig. 3. Illustration of the mechanical (left) and electrophysiological (right) natural and essential boundary conditions.

F ¼ F eFa:

(1)

This multiplicative decomposition inherently results in an incompatible intermediate vector space Ba . In contrast to the total deformation gradient F, neither Fe nor Fa is related to a gradient of any non-linear deformation map (Fig. 2). The active part of the deformation gradient Fa evolves with the transmembrane potential F and carries the information about the underlying actively contracting anisotropic architecture of cardiac tissue through its dependence upon the second-order structural tensors Am, An, and Ak a

b ðF; Am ; An ; A Þ: Fa ¼ F k

(2)

For a contractile material possessing an orthotropic microstructure, the active part of the deformation gradient can be expressed as

 a   a   a  F a ¼ 1 þ lm  1 Am þ ln  1 An þ lk  1 Ak

(3)

a

l a ðFÞ for a ¼ m, k, n denote the potential-dependent where la ¼ b active stretches in the preferred directions. The passive part of the deformation gradient Fe then ensures the compatibility of the total deformation. It is clear from the symmetric form of Fa that the rotations are accounted for by the passive part Fe alone. a

f Nagumo et al., 1962), the current source bI controls the characteristics of the action potential with regard to its shape, duration, restitution, and hyperpolarization along with another variable, the so-called recovery variable r. The evolution of the recovery variable r is described by an additional ordinary differential equation (13)2. Since the recovery variable r governs the local repolarization response of the action potential, we treat it as a local internal variable. Analogous to the momentum balance (4), the excitation equation is subject to the following essential and natural boundary conditions, Fig. 3 (right),

F ¼ F on vS f and h ¼ h on vS h ;

(7)

respectively. Clearly, the surface subdomains vS f and vS h are disjoint, vS f XvS h ¼ B, and complementary, vS ¼ vS f WvS h . The electrical surface flux term h in (7)2 is related to the spatial flux b $n in terms of b through the Cauchy-type formula h : ¼ J 1 q vector q the spatial surface normal n. The material time derivative in the excitation Equation (6) necessitates an initial condition for the potential field at t ¼ t0

F0 ðXÞ ¼ FðX; t0 Þ in B:

(8)

b, s; q Observe that the “hat” symbol used along with the terms b f

and bI indicates the dependency of these variables on the primary fields.

2.3. Constitutive equations 2.2. Governing field equations: strong form The quasi-static mechanical equilibrium is described by the local balance of linear momentum equation in its following spatial form

i h s þB ¼ 0 J div J 1 b

in

B

(4)

s and a given body in terms of the Eulerian Kirchhoff stress tensor b force B per unit reference volume. The operator div[] denotes the divergence with respect to the spatial coordinates x. The dependence of the momentum balance on the primary field variables, the deformation field and the electrical field, is indicated by the s, whose specific form is given in Section 3. Kirchhoff stress tensor b The essential (Dirichlet) and natural (Neumann) boundary conditions, see Fig. 3 (left), 4 ¼ 4

on vS 4

and

t ¼ t

on vS t ;

(5)

conclude the description of the mechanical problem. Obviously, the surface subdomains vS 4 and vS t satisfy the conditions vS ¼ vS 4 WvS t and vS 4 XvS t ¼ B. The spatial surface stress traction vector t on vS t is related to the Cauchy stress tensor through the Cauchy stress theorem t : ¼ J 1 s$n through the outward unit surface normal n, defined on vS.

The solution of the coupled electromechanical problem described by the two field Equations (4) and (6) along with the corresponding boundary and initial conditions (5, 7, 8) requires the knowledge of constitutive equations describing the Kirchhoff stress f b , and the current source bI . Following s, the potential flux q tensor b our recently proposed generalized three-dimensional Hill model (Göktepe et al., 2013a,b), in contrast to the literature (Cherubini et al. (2008); Ambrosi et al. (2011)), we additively decompose the p J and the active part free energy function into the passive part c a c J , (Ask et al., 2012a,b),

p a J ¼ c J ðg; FÞ þ c J ðg; F e Þ:

(9)

The passive part of the free energy function depends solely on the total deformation gradient, while the active part depends both on the deformation and on the potential through the elastic part of the deformation gradient. This additive form leads us to the decoupled stress response

s ¼ bs p ðg; FÞ þ bs a ðg; F e Þ

(10)

where the Kirchhoff stress tensor is obtained by the Doyle-Ericksen formula s: ¼ 2vgJ and the elastic part of the deformation gradient is obtained as Fe ¼ FFa1 from (1). Because of the Eulerian structure


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

5

of the formulation the current metric g is explicitly included in the arguments of the constitutive functions. b is assumed to have the following form The potential flux q

b b ¼ Dðg; q FÞ$Vx F

(11)

in terms of the spatial potential gradient VxF and the deformationb dependent anisotropic spatial conduction tensor Dðg; FÞ that governs the conduction speed of the non-planar depolarization front in three-dimensional anisotropic cardiac tissue. The last constitutive relation that describes the electrical source term of the FitzhugheNagumo-type excitation Equation (6) is assumed to be purely electrical

bI f ¼ bI f ðF; rÞ:

(12)

That is, in this contribution, we exclude mechano-electrical coupling which accounts for the opening of ion channels under the action of deformation (Kohl et al., 1999; Nash and Panfilov, 2004; Göktepe and Kuhl, 2010; Göktepe et al., 2013a,b). This purely electrical form describes the effective current generation due to the inward and outward flow of ions across the cell membrane. The ionic flow is triggered by a perturbation of the resting potential beyond some physical threshold upon the arrival of the depolarization front. Apart from the primary field variables, the recovery variable r, which describes the repolarization response of f the action potential, appears among the arguments of bI in (12). The evolution of the recovery variable r chiefly determines the local shape and duration of the action potential inherent to each cardiac cell and may change throughout the heart. For this reason, evolution of the recovery variable r is commonly modeled by a local ordinary differential equation as in (13)2. From an algorithmic point of view, the local nature of this evolution equation allows us to treat the recovery variable as an internal history variable. This is one of the key features of the proposed formulation that preserves the modular global structure of the field equations (Göktepe and Kuhl, 2009; Göktepe et al., 2010; Göktepe and Kuhl, 2010; Dal et al., 2012, 2013; Wong et al., 2011, 2013). The weak formulation of the continuous equations of coupled cardiac electromechanics (4, 6), their linearization and finite element discretization are left out in this work for the sake of conciseness. For the details of the fully coupled implicit finite element formulation of the problem, we refer to our recent works (Göktepe and Kuhl, 2010; Göktepe et al., 2013a,b).

Fig. 4. The AlievePanfilov model with a ¼ 0.01, g ¼ 0.002, b ¼ 0.15, c ¼ 8, m1 ¼ 0.2, m2 ¼ 0.3. The phase portrait depicts trajectories for distinct initial values 40 and r0 (filled circles) converging to a stable equilibrium point (top). Non-oscillatory normalized time plot of the non-dimensional action potential f and the recovery variable r (bottom).

proposed by Noble (1962), Beeler and Reuter (1977), Luo and Rudy (1991) among many others. The reader is also referred to Keener and Sneyd (1998), Sachse (2004), Pullan et al. (2005), Clayton and Panfilov (2008), Tusscher and Panfilov (2008) for the excellent

3. Specific constitutive model In this section, we set forth the specific constitutive model that is used in the numerical simulation of cardiac dysfunctions presented in Section 4. Specifically, we particularize the equations for f s (10), the the current source bI (12), the Kirchhoff stress tensor b evolution of the active part of the deformation gradient Fa (2), the b (11), and the left and right ventricular pressures potential flux q (plv, prv) as natural boundary conditions (5)2 of the mechanical problem.

3.1. Current source The constitutive modeling of excitable cells can be traced back to the seminal work of Hodgkin and Huxley (1952) on neural cells. Their model, which captures the electric response through four variables, was simplified by Fitzhugh (1961) to a two-variable phenomenological model (Göktepe, 2013). This approach has then been followed by the action potential models of cardiac cells

Fig. 5. Anisotropic architecture of the myocardium. The orthogonal unit vectors f0 and s0 designate the preferred fiber and sheet directions in the undeformed configuration, respectively. The third direction n0 is orthogonal to the latter by its definition n0: ¼ (f0  s0)/jf0  s0j.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

6

review of cardiac electrophysiology models. In this study, we adopt the AlievePanfilov model (Aliev and Panfilov, 1996) that also belongs to the class of two-variable phenomenological models. In their local representation, this class of models involves two ordinary differential equations for the rapidly evolving dimensionless f transmembrane potential vs f ¼ bi ðf; rÞ and the slowly evolving r recovery variable vs r ¼ bi ðf; rÞ. The specific form of these local ordinary differential equations for the AlievePanfilov model is given by f vs f ¼ bi ðf; rÞ ¼ cfðf  aÞð1  fÞ  r f; r ε ðf; rÞ½  r  c fðf  b  1Þ; vs r ¼ bi ðf; rÞ ¼ b

(13)

in terms of the material parameters a, b, and c. The coefficient function b ε ðf; rÞ : ¼ g þ m1 r=ðm2 þ fÞ governs the restitution characteristics of the model through the additional material parameters g, m1, and m2, (Aliev and Panfilov, 1996; Göktepe and Kuhl, 2009). In the context of phenomenological electrophysiology, it is common practice to set up the governing equations and parameters in the dimensionless space. For this purpose, we have introduced the dimensionless transmembrane potential f and the dimensionless time s in (13). The dimensionless quantities are related to their physical counterparts through the following conversion formulas.

F ¼ bf f  df and t ¼ bt s:

(14)

Accordingly, the dimensionless potential f is transformed into the physical transmembrane potential F through the factor b4 and the potential difference d4, which are both in millivolts (mV). Likewise, the dimensionless time s is converted into the physical time t by means of the scaling factor bt, in milliseconds (ms). The f conversion expression between the physical source term bI in (6) f b and the normalized source term i in (13)1 can then be written as

bI ¼ bfbi f : bt f

cardiac cells to the overall pumping function of the heart. It is thus critical that the constitutive equations describing the passive and active stress response of cardiac tissue, as well as the one controlling the conductivity, need to incorporate this inherently anisotropic micro-structure. For the passive response of myocardium (9), we employ the orthotropic hyperelasticity model recently proposed by Holzapfel and Ogden (2009), see also Göktepe et al. (2011) for its threedimensional algorithmic finite element implementation and its extension to quasi-incompressible hyperelasticity. The specific form of this model can be expressed through the following free energy function p ~ ðI ; I ; I ; I Þ c J ðg; FÞ ¼ UðJÞ þ J 1 4m 4n 4k

(16)

where U(J) is the purely volumetric part depending on the volume ~ ðI ; I ; I ; I Þ is defined map J: ¼ det(F) and the orthotropic part J 1 4m 4n 4k as

~ ¼ J ~ ðI ; I ; I ; I Þ J 1 4m 4n 4k ¼

h i o X ai n a1 exp bi ðI4i  1Þ2  1 exp½bðI1  3Þ þ (17) 2b1 2bi i ¼ m;n   i ak h 2 exp bk I4k 1 þ 2bk

in terms of the material parameters a1, b1, am, bm, an, bn, ak, bk and the invariants I1, I4m, I4n, and I4k

I1 : ¼ g : b; I4m : ¼ g : m; I4n : ¼ g : n; I4k : ¼ g : k:

(18)

The Eulerian structural tensors m, n, and k are defined as the push-forward of their Lagrangean counterparts

(15)

The phase diagram in Fig. 4 (top) depicts the solution trajectories of the local ordinary differential equations given in (13) corresponding to distinct initial points 40 and r0 denoted by filled f r circles. The dashed nullclines, where bi ¼ 0 or bi ¼ 0, guide the solution. The diagrams in Fig. 4 (bottom) show the dimensionless potential f and the recovery variable r curves plotted against the dimensionless time s. The action potential is generated by adding external stimulation I ¼ 30 to the right-hand side of (13)1 from s ¼ 30 to s ¼ 30.02. In the algorithmic setting, the recovery variable r, introduced in (13)2, is locally stored as an internal variable at each integration point and we use the backward Euler integration scheme to update the current value of r within a typical time step Ds: ¼ ssn, (Göktepe and Kuhl (2009, 2010)). 3.2. Passive stress response The ventricular myocardium can be conceived as a continuum with a hierarchical anisotropic architecture where uni-directionally aligned muscle fibers are interconnected in the form of sheets. These approximately four-cell-thick sheets that are loosely connected by perimysial collagen can easily slide along each other while being stiffest in the direction of the large coiled perimysial fibers aligned with the long axes of the cardiomyocytes, as depicted in Fig. 5. This heterogeneous but well-organized architecture is of fundamental importance for the successful transduction of the essentially one-dimensional excitationecontraction of individual

m : ¼ FMF T ; n : ¼ FNF T ; k : ¼ FKF T ;

(19)

where the Lagrangean structural tensors

M : ¼ f 0 5f 0 ; N : ¼ s0 5s0 ; K : ¼ symðf 0 5s0 Þ

(20)

reflect the underlying orthotropic micro-structure of the myocardium through the vectors f0 and s0 that denote the preferred fiber and sheet directions of the material micro-structure in the undeformed configuration as depicted in Fig. 5. The explicit form of the s p can be readily obtained through passive Kirchhoff stress tensor b the DoyleeEricksen formula (Doyle and Ericksen, 1956)

b s p ðg; FÞ ¼ J bp g 1 þ ~sp ðg; FÞ

(21)

where the orthotropic part of the passive stress is defined as ~ and b ~sp ðg; FÞ : ¼ 2vg J p : ¼ U 0 ðJÞ ¼ dU=dJ. With the definitions of the invariants (18) at hand, the orthotropic part of the passive Kirchhoff stress tensor can be written as

b s p ðg; FÞ ¼ 2J1 b þ 2J4m m þ 2J4n n þ 2J4k k

(22)

with b: ¼ FG1FT denoting the left Cauchy-Green tensor. Moreover, the deformation-dependent scalar stress coefficients J1, J4m, J4n, and J4k are none other than the scalar derivatives of the passive free energy (17) with respect to the invariants


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

a1 exp½b1 ðI1  3Þ; 2 i h ~ ¼ am ðI  1Þexp bm ðI  1Þ2 ; J4m : ¼ vI4m J 4m 4m i h ~ ¼ an ðI  1Þexp bn ðI  1Þ2 ; J4n : ¼ vI4n J 4n 4n i h ~ ¼ a I exp b I 2 : J4k : ¼ vI4k J k 4k k 4k ~ ¼ J1 : ¼ vI1 J

Table 1 Material parameters of the specific model.

(23)

For the active part of the free energy function (9), we assume the following transversely isotropic function

(24)

e : ¼ g : me in terms of the active modulus h and the invariant I4m e e eT with m : ¼ F MF . This leads to us to the active part of the Kirchhoff stress tensor (10)

 e b s a ðg; F e ; MÞ ¼ 2 h I4m  1 me :

(25)

s a requires the knowledge Calculation of the active stress tensor b of the elastic part of the deformation gradient, which depends on the active part of the deformation gradient. Through the following Ansatz, we assume that Fa implicitly depends on the transmembrane potential through the normalized calcium concentration c  a  l ðcÞ  1 M: Fa ¼ 1 þ b

(26)

We can then express the elastic part of the deformation gradient through Fe ¼ FFa1 that yields the following closed-form expression

 a1  l Ff 0 5f 0 : Fe ¼ F  1  b

(27)

The contraction of myocytes is initiated by calcium released from the sarcoplasmic reticulum into the cytosol. This release is initiated by the excitation of myocytes where the transmembrane potential undergoesa large excursions, see Fig. 4 (bottom). Hence, l is considered to be a function of the normalthe active stretch b ized intracellular calcium concentration c through the following relationship a b l ðcÞ ¼

x 1 þ f ðcÞðx  1Þ

lamax ;

(28)

a

where the parameter lmax controls the maximum contraction. Moreover, the functions f and x are defined as

1 1 þ arctanðb ln cÞ; 2 p f ðc0 Þ  1 x :¼ : a f ðc0 Þ  lmax f ðcÞ : ¼

(29)

Here, cðt0 Þ ¼ c0 ¼ 0 is the initial value of the normalized calcium concentration and b denotes the material parameter controlling the rate of switch. When c ¼ c0 ¼ 0, these functions take a1 the values f ðc0 Þ ¼ 0 and x ¼ lmax . Inserting these initial values a into (28), we obtain l t0 ¼ 1. We describe the intracellular calcium concentration transient during the course of the depolarization and repolarization phases by the following evolution equation

vs c ¼ q f  kc;

Parameter

Description

Equation

a1, b1 am, bm, an, bn, ak, bk a q; k; b; lmax a, b, c g, m1, m2

Isotropic passive stress response Anisotropic passive stress response Active stress response Evolution of active contraction Dynamics of the AP model Restitution properties

(17) (17) (24) (28, 29) (13) (13)

diso, dani

Conduction speed

(32)

q; mat ; mar ; c; pat ; par

Pressure evolution (Signorini)

(33, 34)

Cap, Rc, Rp

Pressure evolution (Windkessel)

(35)

h

3.3. Active stress response and active contraction

a 2 1 e c J ðg; F e ; MÞ ¼ h I4m 1 2

7

(30)

adopted from Pelce et al. (1995). In this constitutive equation, the parameters q and k control the steady-state value of c and its rate of evolution (Göktepe et al., 2013a,b). Specifically, for a given value of the dimensionless action potential f, the steady-state normalized calcium concentration cN is controlled by the ratio q/k. The parameter q governs the initial rate of evolution. This evolution of the normalized calcium concentration c in (30) is integrated locally by using the implicit Euler scheme for the current value of c within a typical time step Ds

c ¼ ðcn þ DsqfÞ=ð1 þ DskÞ:

(31)

3.4. Spatial potential flux b in (11) as We have already introduced the spatial potential flux q a function of the conduction tensor D and the potential gradient VxF. For the model problem, we decompose the second-order conduction tensor D into isotropic and anisotropic parts

D ¼ diso g 1 þ dani m

(32)

in terms of the scalar conduction coefficients diso and dani, where the latter accounts for the faster conduction along the myofiber directions. 3.5. Ventricular pressures The left ventricular pressure plv is modeled by using the Signorini model (Sainte-Marie et al., 2006) for the isovolumetric contraction, isovolumetric relaxation, and filling phases

~Þ plv : ¼ b p0 þ ðb pN  b p 0 Þexp½  expðcq

(33)

~ : ¼ q  q. Here, q denotes the with the normalized flow defined as q flow from the left ventricle and is defined as the negative rate of the left ventricular volume change q : ¼ V_ lv . The parameters q and c control the flow of switch and the sharpness of the switch, respectively. The flow-dependent coefficients b p 0 and b p N are defined as

~Þ : ¼ pat þ mat q ~ b p 0 ðq

and

~Þ : ¼ par þ mar q ~ b p N ðq

(34)

in terms of the atrial pressure pat, the arterial pressure par, and the resistance constants mat and mar. We employ a three-element Windkessel model for the description of the evolution of the left ventricular pressure plv during the ejection phase

p_ lv : ¼

  1 Rc p 1þ q þ Rc q_  lv Cap Rp Cap Rp

(35)


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

8

in terms of the parameters Rc, Rp, and Cap controlling the resistance and compliance characteristics of the hemodynamical system. This evolution equation is numerically integrated through the explicit Euler scheme to obtain the left ventricular pressure plv. Owing to the structural differences between the left and right ventricles, the right ventricular pressure prv is approximated by one-fifth of the left ventricular pressure. The material parameters of the specified model are summarized in Table 1 along with their constitutive descriptions and the equation numbers where they are introduced. 4. Modeling of cardiac dysfunctions This section aims to demonstrate the potential of the proposed approach to model different cardiac dysfunctions through the electromechanical analysis of a generic heart model. For this purpose, we consider infarction, eccentric hypertrophy (ventricular dilation), and concentric hypertrophy (ventricular thickening) as examples for cardiac diseases. In addition to the pathological cases, the healthy heart model is also analyzed. The results obtained from the analyses of the diseased heart models are compared to the normal case. Apart from the left ventricular pressureevolume curves, the cross-sectional snapshots generated at different phases of the cardiac cycle are compared to demonstrate the changes in ventricular function for the different dysfunctions considered in this study. 4.1. Healthy heart In the coupled electromechanical analysis of healthy and diseased cases, a biventricular generic heart model composed of two truncated ellipsoids is utilized. The generic heart model, depicted in Fig. 6 (top), is discretized into 13348 four-node coupled tetrahedral elements connected at 3059 nodes. Fig. 6 (bottom left) illustrates the heterogeneous orientation of the fiber direction f0 varying gradually from 70 in the epicardium to þ70 in the endocardium. The position-dependent orientation of sheets s0 is depicted in Fig. 6 (bottom right). The displacement degrees of freedom on the top base surface are restrained and the whole outer surface of the heart is assumed to be electrically flux-free. The initiation of the electrical excitation is triggered by an elevated initial value of F0 ¼ 10 mV at the upper part of the septum, the wall separating the left and right ventricles. The initial value of transmembrane potential at the remaining nodes are set to the resting potential F0 ¼ 80 mV. For the finite element analysis of the normal heart, we use the values of the material parameters given in Table 2. In the conversion (14, 15), we use the factors b4 ¼ 100 mV, d4 ¼ 80 mV and bt ¼ 12.9 ms. These are calibrated to obtain the physiological action potential response ranging from 80 mV to þ20 mV as suggested in Aliev and Panfilov (1996). For the diseased heart models, presented in Sections 4.2 and 4.3, the values of the selected material parameters are modified to mimic the corresponding dysfunction. Otherwise, the values of the other parameters are assumed to be the same as those used for the normal heart, given in Table 2. 4.2. Myocardial infarction Myocardial infarction is one of the leading causes of heart failure. The main reason of infarction is the imbalance between the demand of myocardium to sustain vital metabolic functions and the blood supply. The limitation of blood flow to a part of myocardium for a long time results in the death of the myocytes due to the limited supply of oxygen and nutrients. If not treated immediately, the infarcted region loses its contractility considerably, becomes

Fig. 6. The geometry and discretization of a generic heart model generated by truncated ellipsoids. All dimensions are in millimeters (top). The position-dependent orientation of the myofibers f0(X) (bottom left) and the sheets s0(X) (bottom right) in B.

stiff, and thus, cannot contribute to the overall pumping function of the heart anymore. The heart tries to restore this loss by making the remaining healthy cardiac muscle cells work harder. This, in turn, results in remodeling and damage of these healthy cells as well. The implications of the infarction-related histological and structural malfunctions in the heart can also be realized through the deviations of the hemodynamic parameters of the pressureevolume curves such as SV, EDPVR, and ESPVR (Fig. 1) from their physiological values. One of the most common clinical symptoms of infarction is a decrease in the end-systolic pressure value during the ejection phase due to the reduced ejection of the blood into the aorta. This leads to more compliant ESPVR and the reduced SV. The latter is a clear indication of a decrease in the amount of the blood volume pumped in each cardiac cycle, thus the reduced cardiac efficiency (Katz (2011)).

Table 2 Values of the material parameters used for the normal heart. Stress response

Active contraction Excitation Conduction Pressure evolution (Signorini)

Pressure evolution (Windkessel)

a1 ¼ 0.496 kPa, b1 ¼ 7.209[] am ¼ 15.193 kPa, bm ¼ 20.417[] an ¼ 3.283 kPa, bn ¼ 11.176[] ak ¼ 0.662 kPa, bk ¼ 9.466[] h ¼ 100 kPa q ¼ 0.1[], k ¼ 0.07[] b ¼ 3½; lamax ¼ 0:4½ a ¼ 0.01[], b ¼ 0.15 [], c ¼ 8 [] g ¼ 0.002 [], m1 ¼ 0.2 [], m2 ¼ 0.3 [] diso ¼ 1 mm2/ms, dani ¼ 0.1 mm2/ms q ¼ 228 mm3 =s; c ¼ 1 ms=mm3 pat ¼ 20 mmHg, par ¼ 80 mmHg mat ¼ 5$103 mmHg ms/mm3 mar ¼ 8$102 mmHg ms/mm3 Cap ¼ 950 mm3/mmHg Rc ¼ 103 mmHg ms/mm3 Rp ¼ 1 mmHg ms/mm3


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

In this study, we consider two models of the infarcted heart with different sizes of infarctions to demonstrate the correlation between the size of infarction and the heart’s electromechanical cardiac activity. In Fig. 7 (bottom), the infarcted regions are shown in red (in the web version), while the blue (in the web version) regions represent the healthy myocardium. In the heart with a smaller infarcted region (30 ), Fig. 7 (bottom left), the infarction extends from q ¼ 120 to q ¼ 150 where the angle q is measured from the x-axis within the xy-plane (Fig. 6). In the heart with a larger infarction (50 ), Fig. 7 (bottom right), the infarcted region is situated between q ¼ 110 and q ¼ 160 . The height of infarction in the both models measures 30 mm. Note that due to the patch recovery technique used to project the material set numbers, the zones corresponding to the boundaries between the healthy and infarcted regions appear in shades of green (in the web version) in Figs. 7e9. As mentioned above, the infarction is comprised of stiffer scar tissue with reduced contractility. To model infarction in our computer models, we assign different values to some material parameters from those used for the healthy cardiac muscle (Table 2). In order to achieve stiffer tissue we magnify the bulk modulus by ten, a the isotropic modulus is altered to a1 ¼ 496 kPa, and we set lmax ¼ 1 to account for the reduced contractility in the infarcted tissue. The influence of infarction and its size on the electromechanical response of the heart is investigated through the electromechanical finite element analyses of the normal heart model and the infarcted heart models with 30 and 50 sizes of infarction, Fig. 7 (bottom). To appreciate the decrease in the contractility of the infarcted cells better, we have taken lateral and longitudinal cross-sectional slices at z ¼ 20 mm and q ¼ 135 , respectively, as shown in Fig. 7 (top). The snapshots of the lateral and longitudinal slices at the systolic and diastolic phases of the cardiac cycle are demonstrated in Figs. 8 and 9, respectively. Here, the first columns contain the snapshots of the normal heart, while the second and third columns demonstrate

9

the snapshots of the infarcted hearts with 30 and 50 sizes of infarction, respectively. In the first columns of Figs. 8 and 9, we observe the contraction and relaxation of the cross-sectional slices in the systolic and diastolic phases of the cardiac cycle for the normal heart. In systole (Fig. 8, rows 1e3) the excitation-induced contraction of anisotropic myocardium results in physiological thickening of the ventricular walls and the overall twist of the heart about the z-axis. Clearly, the electromechanical thickening of ventricular walls leads to the reduction in the ventricular cavity volumes. The physiological wall thickening and the upward motion of the apex during systole can also be appreciated in the longitudinal slices in the rows 1e3 of Fig. 9. Upon repolarization, the thickened walls relax back to their initial shape, the overall heart untwist, and the apex move downward during diastole (Figs. 8 and 9, rows 4e5). As we compare the cross-sectional snapshots generated for the infarcted models in the second and third columns of Figs. 8 and 9 with those plotted for the normal heart, we distinctly notice the reduced contraction and wall thickening in the infarcted regions indicating the decreased overall pumping efficiency of the heart. Moreover, the fairly sharp strain gradients observed at the interface between the healthy and infarcted tissue at end-systole (Figs. 8 and 9, row 3) implies the stress concentrations and excessive deformations in the neighborhood of infarction. Apart from the deformed cross-sectional slices depicted in Figs. 8 and 9, the PV curves of the infarcted hearts are plotted in Fig. 10 in comparison with PV curve of the healthy heart. In the calculation of the left ventricular pressure for the infarcted hearts, the value of the switch flow parameter q in (33) is changed to q ¼ 198 mm3 =ms in the model with 30 infarction and to q ¼ 191 mm3 =ms in the model with 50 infarction to be able to achieve the steep pressure increase in the isovolumetric contraction. Comparing the PV curves in Fig. 10, we observe that the increase in the infarction size leads to a more compliant ESPVR with lower the ventricular pressure in the ejection phase. These changes cause a decrease in the stroke volume, indicating the impaired pumping capability of the heart. It is also worth mentioning that as the infarcted region gets larger, the effect of the infarction on the PV curve becomes more pronounced. 4.3. Concentric and eccentric hypertrophy In this section, we consider cardiac growth as an additional example of cardiac disease besides the infarcted heart models discussed in the preceding section. Cardiac growth can be conceived as remodeling of cardiac tissue in response to chronic overloading associated with the excessive systolic and diastolic wall stresses. On the cell level, maladaptive cardiac growth takes place through sarcomerogenesis, the process of creation and deposition of new sarcomere units (Opie, 2004; Emmanouilides et al., 1994; Kumar et al., 2005). An elevated systolic wall stress (chronic pressure overload) causes sarcomeres to add in parallel. Hence, the cross-sectional area of cardiomyocytes increases. An increased diastolic wall stress (chronic volume overload), however, results in the deposition of sarcomeres in series. This leads to an increase in the length of cardiomyocytes. On the organ level, the former manifests itself in the form of concentric hypertrophy (ventricular wall thickening), while the latter is associated with the eccentric hypertrophy (ventricular dilation).

Fig. 7. The locations of the lateral and longitudinal cross-sectional slices, used to generate the snapshots in Fig. 8 and 9 (top). The positions and dimensions of the infarcted regions (bottom). Two sizes of infarction have been considered. The smaller infarction (30 ) extends from q ¼ 120 to q ¼ 150 and the large one (50 ) is situated between q ¼ 110 and q ¼ 160 . The height of both measures 30 mm.

4.3.1. Generation of hypertrophied heart models In order to investigate the effects of maladaptive cardiac growth on the heart function, we first use the theory of volume growth to generate the concentrically and eccentrically hypertrophied heart models. For this purpose, we adopt the formulation proposed in


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

10

  F e ¼ F  1  1 wg a5a0

(38)

where we have introduced a: ¼ Fa0. Fsssor the stress-producing free energy function, we utilize the orthotropic model (16) in terms of the elastic invariants g c ~ g I e ; I e ; I e ; I e  J ðg; F e Þ ¼ U g ðJ e Þ þ J 1 4m 4n 4k

(39)

where Ug(Je) is the elastic volumetric part depending on the volume ~ g ðI e ; I e ; I e ; I e Þ has the map Je: ¼ det(Fe) and the orthotropic part J 1 4m 4n 4k e ; I e ; I e are calculated as same form as (17) but the invariants I1e ; I4m 4n 4k the appropriate norms (18) of the preferred directions mapped by the elastic part of the deformation gradient Fe (38). The local evolution of the growth variable wg is described by the following ordinary differential equation

  w_ g ¼ b g wg z ;

(40)

Fig. 8. The snapshots of the lateral slices at different phases of the cardiac cycle are generated by the electromechanical finite element analyses of the normal heart, first column, and the infarcted hearts with 30 and 50 sizes of infarction in the second and third columns, respectively. The lateral cross-section is located at z ¼ 20 mm as depicted in Fig. 7 (top).

Göktepe et al. (2010a,b) where the deformation gradient F is multiplicatively decomposed into the elastic deformation gradient Fe and the growth tensor Fg

F ¼ F eFg:

(36)

We consider the following representation of the growth tensor

  F g ¼ 1 þ wg  1 a0 5a0

(37)

in terms of the scalar growth field wg and the generalized unit vector a0 ˛ TxB that coincides with s0 in concentric hypertrophy and with f0 in eccentric hypertrophy, respectively. The rank-one modified growth tensor (37) can be inverted analytically to obtain the elastic part of the deformation gradient

Fig. 9. The snapshots of the longitudinal slices at different phases of the cardiac cycle are generated by the electromechanical finite element analyses of the normal heart, first column, and the infarcted hearts with 30 and 50 sizes of infarction in the second and third columns, respectively. The longitudinal cross-section is positioned at q ¼ 135 as depicted in Fig. 7 (top).


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

11

(dilated) heart models are illustrated in Fig. 11 (bottom) in comparison with the normal heart model (Fig. 11, top). The contour plots of the growth variable wg illustrate the distribution of the concentric and eccentric hypertrophy throughout the generic heart models. In the concentrically grown heart (Fig. 11, bottom left), we observe the significant thickening in the cardiac walls and a slight increase in the left ventricular volume. In the eccentrically grown heart (Fig. 11, bottom right), however, we note a significant increase in the left ventricular volume and almost no change in the thickness of the cardiac walls. These results are qualitatively in accord with the clinical observations reported on the thickened and dilated hearts.

Fig. 10. Left ventricular pressureevolume curves of the generic heart model for healthy and infarcted cases.

where the function b g ðwg Þ ensures the boundedness of growth and governs the speed of the evolution of growth, while the term z acts as a growth criterion; that is, the growth parameter evolves only when z > 0. The specific form of the growth function is defined as



  1 wmax  wg n b g wg ¼ s wmax  1

(41)

in terms of the steady-state growth wmax and the additional parameters s and n governing the rate and non-linearity of the evolution (40). The growth criterion function z plays a key role in modeling different types of cardiac growth stemming from distinct underlying physiological mechanisms. Here, we consider two distinct cases of maladaptive cardiac growth: the pressure overload-driven concentric growth and the volume overloaddriven eccentric growth. These are differentiated from each other through the following distinct functional dependencies of the growth criterion

z ¼ bz ðsÞ ¼ trðsÞ  pcrt ;

  z ¼ ~z le ¼ le  lcrt :

4.3.2. Electromechanical response of grown hearts To investigate the influence of concentric and eccentric hypertrophy on the overall pumping function of the heart, we have conducted fully coupled electromechanical analyses of the thickened and dilated heart models, depicted in Fig. 11 (bottom). In order to account for the reduced contractility accompanying hypertrophy, we reduced the compressibility of the hypertrophied models 33% and kept the other material parameters unchanged as given in Table 2 for the normal heart. In the second and third columns of Fig. 12, we illustrate the contraction and relaxation of the lateral cross-sectional slices in the systolic and diastolic phases of the cardiac cycle in comparison with the slices generated for the normal heart in the first column. The characteristic features such as contraction-induced thickening of the ventricular walls and the overall twist of the heart about its longitudinal axis are commonly observed in the healthy and hypertrophied heart models (Fig. 12). Although the end diastolic volumes (EDV) of the normal heart and the thickened heart models are comparable, we note that the end-systolic volume (ESV) of the thickened heart model is higher than that of the normal heart model, see also Fig. 13. Comparing the end-systolic volumes of the normal and dilated heart models whose EDVs are substantially different, we observe that the ESV volume of the dilated heart is

(42)

While the former accounts for the cardiac wall thickening z ðsÞ, the latter inthrough the stress-dependent criterion z ¼ b corporates the cardiac dilation due to the overstretching of myozðle Þ. In the fibers through the stretch-dependent criterion z ¼ ~ specific growth criteria tr(s) ¼ s:g designates the weighted qffiffiffiffiffi(42), ffi e e refers to the elastic fiber stretch. The pressure and l : ¼ I4f activation of the growth is controlled by the critical weighted pressure pcrt and the critical elastic fiber stretch lcrt in the respective cases. To generate the concentrically grown heart model, we set a0 ¼ s0 in (37) and use the values of the growth parameters s ¼ 0.1 MPa/s, n ¼ 2, wmax ¼ 1.5, pcrt ¼ 0.012 MPa. To produce the eccentrically hypertrophied heart model, we set a0 ¼ f0 in (37) and use the values of the growth parameters s ¼ 0.1 s1, n ¼ 2, wmax ¼ 1.2, lcrt ¼ 1.001. The values of the orthotropic free energy function parameters (39) are taken to be the same as the passive stress response parameters listed in Table 2 for the both models. The left ventricles of the two heart models with the above-specified material parameters are then subjected to a one-cycle pressure transient given in Fig. 6 of Göktepe et al. (2010a), while one-fifth of this pressure is applied to the right ventricles. The concentrically hypertrophied (thickened) and the eccentrically hypertrophied

Fig. 11. The healthy (top), concentrically hypertrophied (bottom left), and eccentrically hypertrophied (bottom right) heart models. The contour plots of the growth variable wg illustrate the distribution of the concentric and eccentric hypertrophy throughout the heart.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

12

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

Fig. 13. Left ventricular pressureevolume curves obtained through the electromechanical analyses of the normal, thickened, and dilated heart models.

some additional parameters related to the contractility of myocytes and the pressure evolution to investigate additional factors leading to the heart failure. 5. Concluding remarks

Fig. 12. The snapshots of the lateral slices at different phases of the cardiac cycle are generated by the electromechanical finite element analyses of the normal heart (first column), the thickened heart (second column), and the dilated heart (third column). The lateral cross-section is located at z ¼ 20 mm as depicted in Fig. 7 (top). The contour plots demonstrate the distribution of the transmembrane potential F throughout the lateral slices.

significantly greater than the ESV of the normal and thickened heart models. This indicates the decreased overall pumping efficiency of the hypertrophied hearts. The qualitative observations made on the ventricular volumes at diastole and systole based on the cross-sectional slices in Fig. 12 can also be made from the pressureevolume curves generated for the normal and hypertrophied heart models in Fig. 13. Despite the substantial changes in the end-diastolic and end-systolic volumes in the dilated heart compared to the normal heart, we note that the end-systolic pressure in the dilated heart model is slightly higher than that of the normal heart. In the thickened heart model, however, we observe a significantly higher end-systolic pressure than that of the normal heart. Our findings obtained form the computational analyses, demonstrated in Figs. 12 and 13, favorably agree with the clinical observations related to the concentric and eccentric hypertrophy (Klabunde, 2005; Katz, 2011). Obviously, one can alter

In this work, we have developed novel computational heart models where the heart diseases related to infarction, eccentric hypertrophy, and concentric hypertrophy have been incorporated. The coupled problem of cardiac electromechanics through the conservation of linear momentum equation and the excitation equation has been formulated in the Eulerian setting. We have modeled the electromechanical coupling with the recently proposed, novel generalized Hill model (Göktepe et al., 2013a,b) where the classical Hill model (Hill, 1938, 1970) has been extended to the three-dimensional finite deformation setting. This approach has allowed us to recover the advantageous features of the active-strain and active-stress-based models within a unified constitutive framework. The proposed coupled approach has been further furnished by the Windkesel-based models to account for the transient pressure evolution within the ventricular chambers during the cardiac cycle. The proposed model has been employed to investigate different pathological cases that include infarction, eccentric hypertrophy, and concentric hypertrophy. The effects of these distinct cardiac dysfunctions on the pressureevolume curves have been investigated. The predictive capacity of our model is illustrated through the coupled finite element simulations of the healthy and diseased heart models. The computational results demonstrate that the simulated PV curves favorably resemble the typical PV curves corresponding to the investigated cardiac dysfunctions, reported in medical literature. Note that the present work focuses primarily on the effect of changes in the contractile and mechanical properties of the infarcted and hypertrophied tissue on the overall efficiency of the heart. Clearly, the conductivity and excitability of diseased tissue may also be affected besides its active contractile and passive compliance properties. Moreover, the duration of the action potential can also be altered in infarcted and hypertrophied tissue. Apart from these local electrophysiological implications, the asymmetrically grown geometry of the heart may cause deviations in the heart electrical axis, and thus the electrocardiograms. In the follow-up works, we plan to incorporate these effects on the


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

electrophysiological properties of the cardiac tissue and investigate their impact on pressureevolume curves and electrocardiograms. Acknowledgments This study was supported by the European Union Seventh Framework Programme (FP7/2007e2013) under Grant PCIG09-GA2011-294161. References Aliev, R.R., Panfilov, A.V., 1996. A simple two-variable model of cardiac excitation. Chaos Solit. Fractals 7, 293e301. Ambrosi, D., Arioli, G., Nobile, F., Quarteroni, A., 2011. Electromechanical coupling in cardiac dynamics: the active strain approach. SIAM J. Appl. Math. 71, 605e621. Ask, A., Menzel, A., Ristinmaa, M., 2012a. Electrostriction in electro-viscoelastic polymers. Mech. Mater. 50, 9e21. Ask, A., Menzel, A., Ristinmaa, M., 2012b. Phenomenological modeling of viscous electrostrictive polymers. Int. J. Non-Linear Mech. 47, 156e165. Beeler, G.W., Reuter, H., 1977. Reconstruction of the action potential of ventricular myocardial fibres. J. Physiol. 268, 177e210. Burkhoff, D., Mirsky, I., Suga, H., 2005. Assessment of systolic and diastolic ventricular properties via pressureevolume analysis: a guide for clinical, translational, and basic researchers. Am. J. Physiol. Heart Circ. Physiol. 289, H501e H512. Cherubini, C., Filippi, S., Nardinocchi, P., Teresi, L., 2008. An electromechanical model of cardiac tissue: constitutive issues and electrophysiological effects. Prog. Biophys. Mol. Biol. 97, 562e573. Clayton, R.H., Panfilov, A.V., 2008. A guide to modelling cardiac electrical activity in anatomically detailed ventricles. Prog. Biophys. Mol. Biol. 96, 19e43. Constantino, J., Hu, Y., Trayanova, N.A., 2012. A computational approach to understanding the cardiac electromechanical activation sequence in the normal and failing heart, with translation to the clinical practice of crt. Prog. Biophys. Mol. Biol. 110, 372e379. Dal, H., Göktepe, S., Kaliske, M., Kuhl, E., 2012. A fully implicit finite element method for bidomain models of cardiac electrophysiology. Comput. Methods Biomech. Biomed. Eng. 15, 645e656. Dal, H., Göktepe, S., Kaliske, M., Kuhl, E., 2013. A fully implicit finite element method for bidomain models of cardiac electromechanics. Comput. Methods Appl. Mech. Eng. 253, 323e336. Doyle, T.C., Ericksen, J.L., 1956. Nonlinear elasticity. In: Dryden, H.L., von Kármán, T. (Eds.), Advances in Applied Mechanics, vol. 4. Academic Press, New York, pp. 53e116. Emmanouilides, G.C., Riemenschneider, R.A., Allen, H.D., Gutgesell, H.P., 1994. Moss and Adams Heart Disease in Infants, Children, and Adolescents. Lippincott Williams & Wilkins. Eriksson, T.S.E., Prassl, A.J., Plank, G., Holzapfel, G.A., 2013a. Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction. Math. Mech. Solids 18, 592e606. Eriksson, T.S.E., Prassl, A.J., Plank, G., Holzapfel, G.A., 2013b. Modeling the dispersion in electromechanically coupled myocardium. Int. J. Numer. Methods Biomed. Eng. 29, 1267e1284. Fitzhugh, R., 1961. Impulses and physiological states in theoretical models of nerve induction. Biophys. J. 1, 455e466. Formaggia, L., Quarteroni, A., Veneziani, A., 2010. Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System. Springer. Go, A.S., Mozaffarian, D., Roger, V.L., Benjamin, E.J., Berry, J.D., Borden, W.B., Bravata, D.M., Dai, S., Ford, E.S., Fox, C.S., Franco, S., Fullerton, H.J., Gillespie, C., Hailpern, S.M., Heit, J.A., Howard, V.J., Huffman, M.D., Kissela, B.M., Kittner, S.J., Lackland, D.T., Lichtman, J.H., Lisabeth, L.D., Magid, D., Marcus, G.M., Marelli, A., Matchar, D.B., McGuire, D.K., Mohler, E.R., Moy, C.S., Mussolino, M.E., Nichol, G., 2013. Executive summary: heart disease and stroke statistics2013 update: a report from the American heart Association. Circulation 127, 143e152. Göktepe, S., 2013. Fitzhugh-nagumo equation, Springer. In: Engquist, B. (Ed.), To appear in Encyclopedia of Applied and Computational Mathematics. Göktepe, S., Abilez, O.J., Kuhl, E., 2010a. A generic approach towards finite growth with examples of athlete’s heart, cardiac dilation, and cardiac wall thickening. J. Mech. Phys. Solids 58, 1661e1680. Göktepe, S., Abilez, O.J., Parker, K.K., Kuhl, E., 2010b. A multiscale model for eccentric and concentric cardiac growth through sarcomerogenesis. J. Theor. Biol. 265, 433e442. Göktepe, S., Acharya, S.N.S., Wong, J., Kuhl, E., 2011. Computational modeling of passive myocardium. Int. J. Numer. Methods Biomed. Eng. 27, 1e12. Göktepe, S., Kuhl, E., 2009. Computational modeling of cardiac electrophysiology: a novel finite element approach. Int. J. Numer. Methods Eng. 79, 156e178. Göktepe, S., Kuhl, E., 2010. Electromechanics of the heart: a unified approach to the strongly coupled excitationecontraction problem. Comput. Mech. 45, 227e243. Göktepe, S., Menzel, A., Kuhl, E., 2013a. The Generalized Hill Model: a Kinematic Approach Towards Active Muscle Contraction (submitted for publication). Göktepe, S., Menzel, A., Kuhl, E., 2013b. Micro-structurally based kinematic approaches to electromechanics of the heart. In: Holzapfel, G.A., Kuhl, E. (Eds.),

13

Computer Models in Biomechanics. From Nano to Macro. Springer ScienceþBusiness Media Dordrecht, pp. 175e187 chapter 13. Göktepe, S., Wong, J., Kuhl, E., 2010. Atrial and ventricular fibrillation: computational simulation of spiral waves in cardiac tissue. Arch. Appl. Mech. 80, 569e 580. Hill, A.V., 1938. The heat of shortening and the dynamic constants of muscle. Proc. R. Soc. Lond. Ser. B Biol. Sci. 126, 136e195. Hill, A.V., 1970. First and Last Experiments in Muscle Mechanics. Cambridge University Press. Hodgkin, A., Huxley, A., 1952. A quantitative description of membrane current and its application to excitation and conduction in nerve. J. Physiol. 117, 500e544. Holzapfel, G.A., Ogden, R.W., 2009. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philos. Trans. Ser. A Math. Phys. Eng. Sci. 367, 3445e3475. Katz, A.M., 2011. Physiology of the Heart. Wolters Kluwer Health/Lippincott Williams & Wilkins Health, Philadelphia, PA. Keener, J.P., Sneyd, J., 1998. Mathematical Physiology. Springer, New York. Keldermann, R., Nash, M., Panfilov, A., 2007. Pacemakers in a ReactioneDiffusion mechanics system. J. Stat. Phys. 128, 375e392. Kerckhoffs, R., Healy, S., Usyk, T., McCulloch, A., 2006. Computational methods for cardiac electromechanics. Proc. IEEE 94, 769e783. Kerckhoffs, R.C., McCulloch, A.D., Omens, J.H., Mulligan, L.J., 2009. Effects of biventricular pacing and scar size in a computational model of the failing heart with left bundle branch block. Med. Image Anal. 13, 362e369. Kerckhoffs, R.C.P., Neal, M.L., Gu, Q., Bassingthwaighte, J.B., Omens, J.H., McCulloch, A.D., 2007. Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 35, 1e18. Kerckhoffs, R.C.P., Omens, J.H., McCulloch, A.D., 2012. Mechanical discoordination increases continuously after the onset of left bundle branch block despite constant electrical dyssynchrony in a computational model of cardiac electromechanics and growth. Europace 14 (Suppl. 5), v65ev72. Klabunde, R.E., 2005. Cardiovascular Physiology Concepts. Lippincott Williams & Wilkins. Klabunde, R.E., 2011. Cardiovascular Physiology Concepts. Lippincott Williams & Wilkins. Kohl, P., Hunter, P., Noble, D., 1999. Stretch-induced changes in heart rate and rhythm: clinical observations, experiments and mathematical models. Prog. Biophys. Mol. Biol. 71, 91e138. Kröner, E., 1960. Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen. Arch. Ration. Mech. Anal. 4, 273e334. Kumar, V., Abbas, A.K., Fausto, N., 2005. Robbins and Cotran Pathologic Basis of Disease. Elsevier Saunders. Lafortune, P., Ars, R., Vázquez, M., Houzeaux, G., 2012. Coupled electromechanical model of the heart: parallel finite element formulation. Int. J. Numer. Methods Biomed. Eng. 28, 72e86. Lee, E.H., 1969. Elasticeplastic deformation at finite strain. ASME J. Appl. Mech. 36, 1e6. Luo, C.H., Rudy, Y., 1991. A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circ. Res. 68, 1501e1526. Nagumo, J., Arimoto, S., Yoshizawa, S., 1962. An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061e2070. Nardinocchi, P., Teresi, L., 2007. On the active response of soft living tissues. J. Elast. 88, 27e39. Nash, M.P., Panfilov, A.V., 2004. Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Prog. Biophys. Mol. Biol. 85, 501e522. Nickerson, D., Nash, M., Nielsen, P., Smith, N., Hunter, P., 2006. Computational multiscale modeling in the IUPS physiome project: modeling cardiac electromechanics-Author bios. Syst. Biol. 50. Niederer, S.A., Smith, N.P., 2008. An improved numerical method for strong coupling of excitation and contraction models in the heart. Prog. Biophys. Mol. Biol. 96, 90e111. Nobile, F., Quarteroni, A., RuizBaier, R., 2012. An active strain electromechanical model for cardiac tissue. Int. J. Numer. Methods Biomed. Eng. 28, 52e71. Noble, D., 1962. A modification of the HodgkineHuxley equations applicable to purkinje fibre action and pacemaker potentials. J. Physiol. 160, 317e352. Opie, L.H., 2004. Heart Physiology: From Cell to Circulation. Lippincott Williams & Wilkins. Panfilov, A.V., Keldermann, R.H., Nash, M.P., 2005. Self-organized pacemakers in a coupled reaction-diffusion-mechanics system. Phys. Rev. Lett. 95, 258104e1e 258014e4. Pelce, P., Sun, J., Langeveld, C., 1995. A simple model for excitationecontraction coupling in the heart. Chaos Solit. Fractals 5, 383e391. Pullan, A.J., Buist, M.L., Cheng, L.K., 2005. Mathematical Modeling the Electrical Activity of the Heart. World Scientific, Singapore. Rantner, L.J., Tice, B.M., Trayanova, N.A., 2013. Terminating ventricular tachyarrhythmias using far-field low-voltage stimuli: mechanisms and delivery protocols. Heart Rhythm Off. J. Heart Rhythm Soc. 10, 1209e1217. Roberts, B.N., Yang, P.C., Behrens, S.B., Moreno, J.D., Clancy, C.E., 2012. Computational approaches to understand cardiac electrophysiology and arrhythmias. Am. J. Physiol. Heart Circ. Physiol. 303, H766eH783. Rossi, S., Ruiz-Baier, R., Pavarino, L.F., Quarteroni, A., 2012. Orthotropic active strain models for the numerical simulation of cardiac biomechanics. Int. J. Numer. Methods Biomed. Eng. 28, 761e788.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021

14

E. Berberog
lu et al. / European Journal of Mechanics A/Solids xxx (2014) 1e14

Sachse, F.B., 2004. Computational Cardiology: Modeling of Anatomy, Electrophysiology, and Mechanics. Springer Verlag. Sainte-Marie, J., Chapelle, D., Cimrman, R., Sorine, M., 2006. Modeling and estimation of the cardiac electromechanical activity. Comput. Struct. 84, 1743e 1759. Stålhand, J., Klarbring, A., Holzapfel, G., 2011. A mechanochemical 3D continuum model for smooth muscle contraction under finite strains. J. Theor. Biol. 268, 120e130.

Tusscher, K.H.W.J.T., Panfilov, A.V., 2008. Modelling of the ventricular conduction system. Prog. Biophys. Mol. Biol. 96, 152e170. Wong, J., Göktepe, S., Kuhl, E., 2011. Computational modeling of electrochemical coupling: a novel finite element approach towards ionic models for cardiac electrophysiology. Comput. Methods Appl. Mech. Eng. 200, 3139e3158. Wong, J., Göktepe, S., Kuhl, E., 2013. Computational modeling of chemo-electromechanical coupling: a novel implicit monolithic finite element approach. Int. J. Numer. Methods Biomed. Eng. 29, 1104e1133.


lu, E., et al., Computational modeling of coupled cardiac electromechanics incorporating cardiac Please cite this article in press as: Berberog dysfunctions, European Journal of Mechanics A/Solids (2014), http://dx.doi.org/10.1016/j.euromechsol.2014.02.021